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Summary  of  Annual  Technical  Report  for 
October  1, 1982  -  September  30,  1983 

The  objectives  of  the  research  being  conducted  under  the  current  2  year 
contract  are  to:  (1)  Develop  and  test  methods  of  discrimination  in  the  regional 
and  teleseismic  distance  range  using  physical  source  parameter  discriminants, 
(2)  Pursue  theoretical  and  observational  studies  of  seismic  sources;  (3)  To 
develop  methods  of  theoretical  seismogram  synthesis  in  the  near,  regional  and 
teleseismic  distance  ranges  for  structure  and  source  definition;  (4)  To  develop 
and  apply  advanced  signal  processing/analysis  methods  for  discrimination  and 
explosion  yield  estimation  studies  and;  (5)  to  pursue  near  field  studies  of  explo¬ 
sions  and  earthquakes  for  detailed  source  definition. 

In  this  report  ire  describe  specific  research  results  in:  (l)  The  theoretical 
basis  for  seismic  event  discrimination  and  for  the  inference  of  physical  parame¬ 
ters  describing  seismic  sources;  (2)  Observations  of  body  wave  magnitude  resi¬ 
duals  in  the  U.S.  and  Southern  Canada  and  the  implications  for  yield  estimation 
at  non-U.S.  and/or  U.S.  test  sites;  (3)  Inferences  of  magnitude-yield  relation¬ 
ships  for  both  U.S.  and  U.S.S.R.  test  areas,  with  applications  of  these  relations 
to  U.S.  and  Russian  nuclear  test  observations,  to  obtain  new  yield  estimates:  (4) 
Studies  of  a  variety  of  physical  situations  which  would  produce  seismic  observa¬ 
tions  for  an  earthquake  making  it  appear  "anomalous"  (i.e.,  explosion-like) 
when  viewed  from  the  perspective  of  "standard"  discrimination  methods,  and 
similarly  physical  situations  were  described  for  which  explosions  could  appear 
earthquake-like,  and  which  could  also  lead  to  large  yield  estimate  errors. 

Some  of  the  major  conclusions  drawn  from  the  results  described  are:  (l) 
That  earthquake-explosion  discrimination  by  m*  —  Ms,  at  large  to  moderate 
magnitude  levels,  is  due,  in  part,  to  bas  ic  Mftf^&rFcfes  in  't1hesfif^FW3(5J^hcy 
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radiation  from  the  two  source  types,  as  predicted,with  explosions  producing 
more  high  frequency  seismic  radiation  than  earthquakes  with  the  same  low  fre¬ 
quency  radiation.  Further,  because  of  this  same  spectral  difference  between 
the  two  types  of  seismic  events,  spectral  discrimination,  using  variable  fre¬ 
quency  magnitude  (VFM)  methods,  is  very  effective  in  distinguishing  and  identi¬ 
fying  the  two  event  types  in  the  regional  distance  range  and  can  be  applied  to 
very  low  magnitude  levels  (below  m.*  -  3).  (2)  Deviations  of  observed  single  sta¬ 
tion  body  wave  magnitudes  from  the  mean  observed  magnitude  can  be  very 
Large,  of  the  order  of  one  magnitude  unit  from  the  mean,  and  these  deviations 
can  also  be  strongly  dependent  on  the  azimuth  from  the  station  to  the  source, 
with  variations  as  large  as  .5  magnitude  units  occurring  as  a  function  of 
azimuth.  These  variations  are  considered  to  be  due  to  both  focusing  (defocus- 
ing)  effects  and  variations  of  anelastic  absorption  between  points  of  observa¬ 
tion.  The  azimuthal  variation  at  0B2NV,  at  the  Nevada  test  site,  is  .22  magni¬ 
tude  units,  while  at  RKON,  on  the  Canadian  shield,  the  variation  with  azimuth  is 
.45  units;  suggesting  that  use  of  these  sites  as  analogs  for  the  differences 
between  NTS  and  Russian  test  sites  requires  additional  assumptions  about  how 
the  strong  azimuthal  variations  are  to  be  interpreted.  If,  however,  azimuthal 
variations  are  averaged  and  these  stations  are  used,  then  the  predicted  differ¬ 
ence  in  m*  observations  from  Russian  tests  and  NTS  tests  should  be  .31  magni¬ 
tude  units,  with  NTS  having  the  lower  average  magnitude  if  the  observing  sta¬ 
tions  provide  a  reasonably  uniform  azimuth  sampling  of  both  test  sits.  (3)  Rus¬ 
sian  test  site  body  wave  magnitude-yield  relations,  adjusted  from  NTS  derived 
magnitude-yield  relations  to  account  for  a  transmission  bias  of  slightly  more 
than  .3  magnitude  units  between  NTS  and  Russian  sites,  give  estimated  yields 
for  Russian  tests  (from  1968  to  1977)  that  are  reasonably  consistent  with  yields 
estimated  using  "universal"  surface  wave  magnitude-yield  relations.  Further, 
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ihe  estimated  yields  for  Russian  tests  appear  to  be  within  the  limitations  of  the 
test  ban  agreement.  That  is,  within  the  uncertainties  of  the  estimation 
methods,  all  the  yield  estimates  obtained  are  less  than,  or  about  equal  to,  150 
lcT.  (4)  Seismically  anomalous  events  can  (and  do)  occur  under  a  variety  of 
conditions,  but  it  is  possible  to  properly  identify  these  events  when  appropriate 
seismic  procedures  are  employed,  with  these  methods  being  spectral 
discrimination  methods  which  can  be  implemented  for  routine  analysis. 


I.  Introduction 

The  objectives  of  the  research  being  conducted  under  the  current  2  year 
contract  are  to:  (l)  Develop  and  test  methods  of  discrimination  in  the  regional 
and  teleseismic  distance  range  using  physical  source  parameter  discriminants, 
(2)  Pursue  theoretical  and  observational  studies  of  seismic  sources;  (3)  To 
develop  methods  of  theoretical  seismogram  synthesis  in  the  near,  regional  and 
teleseismic  distance  ranges  for  structure  and  source  definition;  (4)  To  develop 
and  apply  advanced  signal  processing/analysis  methods  for  discrimination  and 
explosion  yield  estimation  studies  and;  (5)  to  pursue  near  field  studies  of  explo¬ 
sions  and  earthquakes  for  detailed  source  definition. 

In  this  report  we  describe  some  of  our  work,  and  results,  on  the  theory  of 
earthquake  source  representation,  as  it  relates  to  earthquake-explosion 
discrimination.  In  this  regard,  Section  II  and  the  Appendix  1  provide  a  descrip¬ 
tion  of  a  major  part  of  our  comprehensive  study  of  the 
theoretical/observational  basis  for  earthquake-explosion  discrimination.  In 
addition,  in  Section  III,  we  summarize  our  observational  results  pertaining  to 
measurements  of  magnitude  variations  for  the  U.S.  and  Southern  Canada,  and 
discus  the  relevance  of  these  observations  to  magnitude  bias  due  to  upper  man¬ 
tle  attenuation,  as  well  as  the  relationship  of  these  results  to  yield  estimation 
for  USSR  underground  tests.  In  a  related  study,  described  in  Section  IV,  we 


•*_  *’  *  •  *,•  *•’■**  *  •  *  •  *  •  '  •  *  •  *  »  m ''  ,  *  ,  «  ,*•  ,  •  \  *.  *,  «  *  ‘  *  -  •  "  k  *  *  .  1 


address  the  effects  of  variable  coupling  and  attenuation  as  they  relate  to  yield 
estimation  for  underground  nuclear  tests  within  the  U.S.  and  the  U.S.S.R.,  and 
then  apply  our  inferences  to  the  Russian  and  U.S.  test  data  in  order  to  obtain 
yield  estimation  results.  In  Section  V  we  describe  some  of  our  investigations  of 
the  conditions  under  which  "anomalous"  events  will  occur,  specifically  for 
-Ms  discrimination  and  VFM  (variable  frequency  magnitude)  discrimination. 

K.  Discrimination  of  Seismic  Events 

Several  research  projects  have  been  pursued  during  the  1982-83  contract 
period.  In  this  report  we  will  emphasize  our  work  on  a  comprehensive  study  of 
the  theoretical  foundations  of  seismic  event  discrimination  (Archambeau  and 
Evernden,  1984).  These  results  are  the  first  part  of  a  comprehensive  study  of 
discrimination,  with  the  complete  study  involving  an  integration  of  previous 
work  on  source  theory,  wave  propagation  theory  and  earth  structure,  including 
new  material  and  results  obtained  under  current  AFOSR  support.  This  work  is 
an  attempt  to  provide  a  complete  explanation  of  empirically  substantiated 
discrimination  methods,  and.  in  addition,  to  provide  a  basis  for  less  verified  or 
new  discrimination  methods.  In  the  complete  study,  the  theoretical  predictions 
are  compared  to  relatively  large  sets  of  observed  data,  in  order  to  substantiate 
the  predictions  and  confirm  the  models.  Further,  inferences  of  network  capa¬ 
bility  for  the  detection  and  identification  of  explosions,  including  decoupled 
explosions,  are  also  considered  in  an  attempt  to  provide  an  accurate  and  realis¬ 
tic  evaluation  of  U.S.  capability  for  detection  and  identification  of  Eurasian 
events  under  a  variety  of  conditions.  Finally,  consideration  is  also  given  to 
anomalous  events,  particularly  anomalous  explosions  and  the  conditions  under 
which  it  is  expected  that  such  events  can  be  generated.  In  this  regard,  joint 
discrimination  procedures  are  described  that  can  be  used  to  identify  such 


The  first  part  of  this  study  is  included  in  the  Appendix  1.  wherein  we  pro¬ 
vide  the  theoretical  foundations  for  quantitative  descriptions  of  earthquake 
generated  seismic  fields,  along  with  specific  spectral  predictions,  scaling  laws, 
etc.  These  results  form  part  of  the  framework  within  which  we  interpret 
discrimination  methods,  the  other  part  being  a  similar  quantitative  representa¬ 
tion  of  explosion  generated  seismic  fields.  Through  the  use  of  these  theoretical 
predictions  we  are  then  able  to  interpret  and  understand  observed  discrimina¬ 
tory  differences  between  earthquakes  and  explosions  (Evernden  and  Archam- 
beau,  1984).  (These  detailed  comparisons  of  theory  and  observation  will  be  dis¬ 
cussed  in  our  next  report.) 

We  are  also  concerned  with  seismic  source  discrimination  through  the  use 
of  inferred  physical  parameters  for  seismic  sources.  In  this  investigation  we 
consider  it  necessary  to  use  formal  inversion  methods  to  obtain  estimates  of 
source  parameters,  such  as  stress  drop,  rupture  length  and  rupture  velocity. 
In  this  regard,  inversion  methods  provide  the  necessary  automated  and  objec¬ 
tive  methods  that  must  be  used,  and  we  have  therefore  investigated  approaches 
which  involve  moment  tensors  as  well  as  new  approaches  to  source  inversion 
based  on  vector  multipole  formulations  of  the  seismic  radiation  fields  from 
explosions  and  earthquakes  (Archambeau  and  Scales,  1984). 

In  this  work  it  has  been  shown  that  the  conceptual  and  physical  basis  for 
Gilbert’s  moment  tensor  representation  (Gilbert,  1971)  is  faulty  and  that  its  cri¬ 
ticism  by  Backus  and  Mulcahy  (I97d)  is  also  in  error.  Further,  it  is  found  that 
the  stress  glut  formulation  and  concept  is  without  physical  or  mathematical 
(t.e.,  logical)  foundation.  An  alternative  approach  for  the  inversion  of  seismic 
data  for  source  properties  is  given  and  it  is  shown  that  a  vector  multipole 
expansion  of  the  seismic  radiation  field  is  appropriate,  in  that  it  provides  a 
desired  linear  inversion  for  these  coefficients  from  the  observed  data,  as  do 


moment  tensor  components.  In  addition,  the  physical  parameters  defining  the 
source,  such  as  spatial  stress  changes,  failure  rate  and  failure  zone  dimensions, 
are  explicitly  related  to  the  vector  multipole  coefficients  and  it  is  shown  how 
these  parameters  can  be  obtained  from  the  multipole  coefficients  through  a 
non-linear  inversion  procedure. 

In  addition  to  a  formal  inversion  approach,  we  have  previously  used  relaxa¬ 
tion  source  models  in  earth  structure  models  to  provide  a  basis  of  interpreta¬ 
tion  of  the  »n.j  and  Ms  data  from  earthquakes,  in  terms  of  failure  zone  dimen¬ 
sions  and  stress  drop  (Archambeau  et  al.,  19B3).  In  particular,  we  studied 
recorded  earthquakes  in  the  Aleutian  region.  During  the  present  contract 
period  this  work  was  finalized  and  is  being  submitted  for  publication.  We  have 
begun  a  similar  study  of  the  magnitude  data  for  the  Japan-Kamchatka  region 
during  the  current  contract  period  and  expect  to  finalize  this  work  by  the  end 
of  the  contract  period.  Finally,  we  have  begun  studies  of  small  earthquakes  and 
explosions  from  the  California-Nevada  region  using  spectral  methods  to  infer 
source  stress-drop,  dimension,  etc.  (Wilson  and  Archambeau,  1984).  One  of  the 
principal  results  is  that  the  average  stress  drops  observed  for  large  numbers  of 
individual  events  show  a  definite  dependence  on  rupture  length,  with  small 
events  giving  much  higher  average  stress  drops  than  larger  (longer)  ones.  This 
result  strongly  suggests  that  many  smaller  events  involve  failure  at  highly 
stressed  single  asperities,  without  failure  zone  propagation  continuing  very  far 
outside  the  single  asperity  zone.  On  the  other  hand,  larger  events  appear  to 
involve  failure  connecting  several  or  a  large  number  of  highly  stressed  asperi¬ 
ties  separated  by  zones  of  low  tectonic  stress.  In  such  a  picture,  explosions 
appear  as  small  rupture  dimension  events  with  very  high  stress  drops  (above 
the  4-5  kilobar  level),  while  earthquakes  of  comparable  effective  failure  zone 
lengF  s  have  ess  drops  significantly  below  the  low  kilobar  level  (below  1  kilo- 


bar)  and  generally,  as  a  population,  have  effective  rupture  dimensions  that  are 
much  larger  than  explosions.  Hence,  a  "physical  parameter  space"  discrim¬ 
inant,  involving  stress  drop  vs  maximum  failure  zone  dimension,  serves  as 
another  approach  to  event  identification. 

During  the  FY  19B2-B3  contract  period  we  also  continued  our  studies  of 
spectral  discriminants,  in  particular,  the  variable  frequency  magnitude  discrim¬ 
inant  employed  in  the  study  of  Eurasian  event  data  set.  This  latter  effort  was 
jointly  conducted  with  Systems,  Science  and  Software.  The  current  effort 
involves  finalization  of  this  work,  along  with  applications  of  new  modeling  tech¬ 
niques,  of  source  and  propagation  effects,  applied  to  provide  a  more  detailed 
understanding  of  the  observational  results.  The  essential  result  obtained,  how¬ 
ever,  is  that  VFM  discrimination  is  very  effective  down  to  very  low  magnitude 
levels  (at  least  to  mb  -  3.0)  and  that  its  effectiveness  is  due  to  the  fact  that 
most  explosions  with  the  same  low  frequency  spectral  level  as  an  earthquake 
will  have  a  higher  "corner  frequency"  than  the  earthquake,  as  well  as  a  lesser 
P-wave  spectral  decay  rate  at  frequencies  above  the  corner  frequency.  Thus, 
this  results  in  substantially  greater  high  frequency  generation  for  an  explosion 
compared  to  an  earthquake  with  the  same  low  frequency  spectral  level  for  the 
direct  P-wave,  and  this  difference  is  exploited  by  the  VFM  method  to  provide  the 
observed  discrimination  capability.  The  method  is  shown  to  be  very  effective  in 
regional  distance  ranges  where  relatively  high  Q  paths  are  available,  particu¬ 
larly  when  very  high  frequency  magnitudes  (up  to  10  Hz  and  beyond)  are  used. 

III.  Magnitude  Bias  Estimates  for  Yield  Estimation 

This  work  (Butler,  1983)  involves  body  wave  magnitude  measurements  at  1 
Hz  from  a  large  set  of  selected  events  at  all  available  U.S.  stations,  as  well  as  at 
some  stations  in  Canada.  At  a  number  of  stations,  including  0B2NV  at  the 
Nevada  test  site,  magnitude  data  was  obtained  at  3  azimuths,  with  each  azimuth 


data  set  analyzed  separately.  The  results  of  this  analysis  are  summarized  in 
Figures  1  and  2.  The  conclusions  that  are  apparent  from  the  results  of  this 
study  are  that  the  relative  body  wave  magnitude  variations  in  the  U.S.  are  large, 
with  the  range  running  at  least  from  +.4  to  -.4  magnitude  units.  Further,  the 
spatial  variation  is,  in  many  areas,  extremely  rapid,  suggesting  focusing  and 
defocusing  effects.  In  this  regard,  the  variations  of  6mb  with  azimuth  at  a  par¬ 
ticular  station  can  be  very  large,  with  variations  of  as  much  as  .4  magnitude 
units  occurring,  and  with  no  apparent  correlation  with  tectonic  province  for 
the  large  variations.  In  particular,  the  Nevada  test  site  station  0B2NV  (on  a 
granite  stock),  shows  an  azimuthal  variation  in  6mb  of  .22  magnitude  units,  and 
the  site  at  RKON,  which  has  been  used  as  a  station  site  comparable  to  Eurasian 
nuclear  test  sites,  shows  an  azimuth  variation  of  .45  magnitude  units,  with 
larger  magnitudes  to  the  north  than  to  the  south  and  northwest.  These  partic¬ 
ular  results  strongly  suggest  that  magnitude  bias  estimates,  using  0B2NV  and 
RKON  as  representative  NTS  and  Russian  test  site  analogs,  may  be  misleading 
and  at  the  least  ambiguous  since  different  azimuthal  sampling  will  give  very  dif¬ 
ferent  results,  [e  g.,  One  could  get  values  for  the  difference,  6mb  (RKON)  -  6mb 
(0B2),  anywhere  in  the  range  from  .74  to  .07  depending  on  azimuth  sampling.)  If 
an  average  of  the  three  azimuths  is  used,  then  6mb  (0B2)  -  .14,  while  6mb 
(RKON)  -  .45,  so  that  5mb  (RKON)  -  6mb  (0B2)  -  .31.  However,  blind  use  of  such 
an  average  in  estimating  magnitude  bias  for  test  sites  is  certainly  suspect, 
unless  one  can  show  that  the  averaging  removes  local  focusing-defocusing  vari¬ 
ations  completely  and  that  the  average  represents  the  basic  anelastic  absorp¬ 
tion  effect  plus  any  regional  focusing-defocusing  due  to  deep  (crust-upper 
mantle)  structure.  Even  so,  one  would  expect  fairly  wide  ranging  fluctuations 


Figure  1.  Distribution  of  stations  and  relative  amplitude  (<5mt,)  observations  for  the  U.S.  and  Southern  Canada 
At  some  of  the  stations  the  data  could  be  analyzed  with  respect  to  three  distinct  azimuths ,  the 
results  of  which  are  indicated  by  three  values  arranged  at  appropriate  azimuths  for  these  stations. 
(Butler,  1983) 


Figure  2.  Contours  of  relative  body  wave  magnitude  for  the  U.S.  and  So.  Canada.  Station  locations  providing  the 

data  are  Indicated.  Azimuthal  data  (3  directions)  was  available  at  stations  indicated  by  the  solid  circles. 
The  azimuthal  variations  observed  at  these  stations  was  taken  into  account  in  the  contouring.  Where  station 
coverage  is  sparse  the  contours  are  questionable,  and  must  be  viewed  more  as  predictions  rather  than 


Basin  and  Range  Province,  and  its  anomalously  high  values  relative  to 

nearby  stations  probably  represents  strong  focusing  associated  with  the  gran¬ 
ite  intrinsive  at  the  site. 

These  results  imply  that  extremely  accurate  yield  estimates  (to  within,  say, 
a  10%  uncertainty)  will  most  likely  require  detailed  site  information  in  view  of 
the  rapid,  and  somewhat  unexpected,  variations  present  within  geologic  pro¬ 
vinces  and  sub-provinces.  However,  correlations  of  these  observations  with 
more  readily  available  geophysical  data  ( e.g .,  gravity,  heat  flow,  P  delays,  cru¬ 
stal  thickness,  near  surface  geology,  Pn  velocities,,  elevation,  etc.)  may  offer 
some  hope  of  deducing  magnitude  variations  for  "random"  sites. 

IV.  Teleseismic  Yield  Estimation  with  Applications  to  U.S.  and  U.S.S.R  Test  Data 

The  necessity  of  obtaining  accurate  and  rather  precise  yield  estimates 
from  seismic  data  places  a  rather  extreme  burden  on  our  methods  of  measure¬ 
ment  and  data  analysis  and,  as  well,  on  the  basis  of  interpretation,  of  such  data. 
As  is  well  known,  a  large  number  of  near  site  environmental  factors,  such  as 
detailed  material  properties  at  the  explosion  site  and  the  nature  of  tectonic 
stress  fields,  are  capable  of  strongly  perturbing  the  seismic  field,  and  so  this 
state  of  affairs  makes  accurate  estimates  of  any  single  parameter  such  as  yield 
extremely  difficult,  if  not  impossible,  without  a  detailed  understanding  of  all 
these  (other)  effects  and  the  means  of  taking  them  into  account  in  the  analysis 
of  the  observed  data.  In  addition,  we  are  also  well  aware  that  variations  in  crust 
and  upper  mantle  properties,  including  both  structure  and  absorption  proper¬ 
ties,  can  lead  to  very  strong  fluctuations  in  the  observed  seismic  data,  resulting 
in  biasing  between  different  test  sites.  Until  we  have  reasonably  accurate 
determinations  of  these  characteristics  for  potential  underground  test  areas,  it 
will  be  difficult  to  obtain  yield  estimates  that  are  reliable  to  within  a  factor  of 
two.  Furthermore,  the  effects  of  anelastic  absorption  and  scattering  (focusing 


and  defocusing)  are  certainly  frequency  dependent  and  without  fairly  precise 
knowledge  of  the  form  of  the  frequency  dependence  for  each,  and  the  relative 
magnitudes  of  these  effects  at  particular  test  locations,  it  is  unlikely  that  very 
precise  yield  estimates  can  be  obtained  because  of  this  uncertainty  alone. 
Finally,  it  seems  clear  that  much  more  precisely  defined  seismic  magnitudes 
are  required  (or  perhaps  a  different  measurement  or  group  of  measurements 
which  are  carefully  and  precisely  defined)  in  order  to  reduce  observational  sub¬ 
jectivity  and  associated  variations  in  the  measured  data. 

In  addition  to  the  effects  of  anelastic  attenuation,  it  is  clear  that  variations 
in  material  properties  produce  changes  in  the  explosion  coupling  and  so  will 
produce  different  seismic  wave  fields  for  the  same  explosive  yield.  Figures  3 
and  4,  from  Archambeau  and  Evernden  (1983),  show  body  and  surface  wave 
data  (from  Marshall  et  al.)  which  illustrates  the  approach  required  to  obtain 
meaningful  yield  estimates.  In  particular,  we  have  used  high,  intermediate  and 
low  coupling  magnitude-yield  curves  to  interpret  this  data,  with  these  curves 
being  based  on  results  derived  from  explosion  modeling  in  different  types  of 
materials.  Further,  the  mb  versus  yield  curves  shown  are  segmented  into  sec¬ 
tions  of  different  constant  slope,  this  change  of  slope  occurring  at  large  yields 
when  the  corner  frequency  for  the  explosive  generated  P  wave  is  less  than  1  Hz. 
The  upper  curve,  for  non-NTS  explosions  in  high  coupling  media,  applies  to  Rus¬ 
sian  event  data  as  well  as  U.S.  Aleutian  explosions.  This  curve,  and  associated 
data,  indicates  a  clear  shift  above  the  high  coupling  NTS  curve,  of  about  .25 
magnitude  units.  This  then  is  a  bias  effect,  which  arises  because  of  the  differ¬ 
ences  between  the  NTS  site  transmission  characteristics  and  the  less 
attenuated  transmission  from  other  sites. 

The  data  show  considerable  scatter  about  the  magnitude  yield  curves,  and 
this  is  considered  to  be  due  to:  deviations  from  the  mean  material  properties 
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3.  Observed  body  wave  magnitude  versus  yield  data  (from  Marshall  et  al.) 
interpreted  using  three  coupling  level  sq,  vs  Y  curves,  which  are  based 
on  modeling  results.  Note  that  yield  curves  generally  consist  of  two 
segments  of  differing  constant  slope,  these  being  due  to  magnitudes 
measured  above  the  corner  frequency  for  large  yield  events,  and  below 
the  corner  frequency  for  intermediate  and  low  yield  explosions. 


Me  vorsus  yield 

*  ( Marshall,  Springer  ft  Rodson,  1979) 

•  *  Announced  yields  (NTS) 

G  •  Announced  yields  (Non  NTS) 

X  •  Announced  yields,  USSR  PNE 
’  T  ■  Tuff 
R  •  Rhyohite 
6  •  Granite 
A  ■  Alluvium 

S  *  Salt 
So  *  Sandstone 
.  C  *  Cloy 


YIELD  (kt) 

Figure  A.  Observed  surf see  wsve  magnitude  versus  yield  data  for  NTS,  non-NTS  and 
USSR  events  with  announced  yields.  The  data  is  interpreted,  as  was  the 
■b  VB  Y  data,  using  three  coupling  level  Mg  vs  Y  relations  which  are 
based  on  modeling  results.  All  world  wide  data  appears  to  be  quite  fit 
by  these  three  coupling  level  curves  of  constant  slope. 


(implicitly)  assumed  for  each  of  the  three  coupling  level  curves;  local  variations 
in  the  transmission  characteristics  arising  from  the  shifts  in  explosion  loca¬ 
tions;  and,  finally,  variations  in  mb  and  Ms  values  due  to  tectonic  release 
effects.  The  latter  effect  has  been  interpreted  as  being  particularly  strong  for 
the  Rulison  and  Rio  Blanco  events  (Murphy,  Archambeau  and  Shah,  1983),  with 
the  values  being  considerably  reduced  by  tectonic  release  having  a  double 
couple  or  fault  plane  equivalent  corresponding  to  a  45°  plane  with  normal 
equivalent  faulting  at  the  explosion  hypocenter.  (The  actual  tectonic  release 
mechanism  is  thought  to  be  stress  relaxation  in  a  prestressed  zone  around  the 
(nearly  spherical)  shatter  zone  created  by  the  explosion  itself.) 

As  examples  of  the  application  of  these  current,  somewhat  preliminary 
results,  to  Russian  explosion  test  site  data,  Figure  5  shows  the  inferred  magni¬ 
tude  yield  curves  for  Russian  (shield)  explosions  with  announced  yields,  along 
with  French  and  U.S.  data  from  the  Aleutians,  the  latter  two  data  sets  con¬ 
sidered  to  be  comparable  to  data  from  Russian  test  sites  in  terms  of  roughly 
equivalent  transmission  characteristics.  In  particular,  the  anelastic  absorption 
and  other  transmission  characteristics  are  considered  to  be  comparable.  The 
empirical  curves  shown  on  this  plot  are  obtained  from  the  NTS  curves  shown 
earlier,  with  a  .3  increase  in  the  mb  value  for  a  given  yield.  This  value  is  based 
on  the  average  5mb  value  between  an  NTS  test  site  (0B2NV)  and  the  shield-like 
station  at  Red  Lake  Ontario  (RKON),  the  result  being  obtained  using  the  data 
shown  in  Figures  1  and  2.  this  value  is  also  consistent  with  the  rough  magni¬ 
tude  of  the  inferred  body  wave  magnitude  difference  between  NTS  and  Aleutian 
events  (.25  mb  units)  indicated  in  Figure  3.  The  fit  of  the  empirical  curves 
shown  in  Figure  5  to  the  observed  data  is  considered  to  be  very  good,  with  the 
high  and  intermediate  coupling  curves  being  in  agreement  with  the  observed 
data  in  a  sensible  and  physically  understandable  manner. 
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Figure  5  .  U.S.  (NTS)  body  wave  aagnltude  versus  yield  curves  adjusted  for 

differences  in  anelastic  absorption  between  NTS  and  USSR  test  sites, 
coopered  to  observed  USSR  events  having  announced  yields.  The  lower  set 
of  three  curves  are  the  universal  Mg  vs  Y  curves  with  non -NTS  surface 
wave  aagnltude  data.  Aleutian  and  French  test  data  are  also  shown  for 
cooparlson.  The  NTS  body  wave  curves  for  high,  intermediate  and  low 
coupling  have  all  been  shifted  upward  by  .3  aagnltude  units. 

(Data  fron  Marshall  et  al.) 
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The  consequences  of  yield  estimation  based  on  magnitude-yield  curves 
incorporating  sensible  corrections  for  strong  biasing  effects,  such  as  differen¬ 
tial  absorption,  variable  coupling  and  tectonic  release  are  indicated  in  Figures 
6  and  7.  Thus,  with  an  uncorrected  application  of  the  NTS  derived  magnitude- 
yield  relations  to  Russian  test  data  one  obtains  yield  estimates  from  conven¬ 
tional  mb  and  Ms  data  as  indicated  in  Figure  6.  As  noted  in  the  figure  caption, 
there  are  several  reasons  why  one  expects  such  estimates  to  be  systematically 
in  error.  On  the  other  hand,  Figure  7  shows  estimates  for  the  same  set  of  Rus¬ 
sian  explosions  when  a  bias  correction  in  mb  values  is  applied  to  the  NTS  magni¬ 
tude  curves  so  as  to  account  for  the  different  transmission  losses  that  are 
inferred  between  NTS  and  the  Russian  sites.  Clearly,  quite  a  different  picture  of 
Russian  test  activity  is  implied  by  Figure  7,  compared  to  that  inferred  from  Fig¬ 
ure  6. 

We  expect  that  the  residual  scatter,  about  the  one  to  one  slope  line 
corresponding  to  yield  values  which  are  the  same  from  both  mb  and  Ms  esti¬ 
mates.  is  mainly  due  to  the  unaccounted  effects  of  tectonic  release  and  local 
variations  of  transmission  characteristics^which  fluctuate  about  the  mean  for 
the  Russian  sites.  In  our  future  work  we  intend  to  attempt  to  significantly 
improve  upon  these  yield  estimates  by  more  fully  and  accurately  accounting  for 
all  three  effects,  and  to  provide  a  strong  justification  for  the  "corrections"  that 
are  used. 

The  yield  estimates  for  U.S.  underground  tests,  before  and  after  the  limited 
test  ban  agreement,  were  obtained  using  the  NTS  magnitude-yield  curves  shown 
earlier  and  these  estimates  may  be  compared  to  the  results  for  the  Russian 
tests.  In  this  regard,  the  U.S.  test  data  gives  yield  estimates  with  a  distribution 
similar  to  the  Russian  results  for  the  1977-81  period.  The  somewhat  greater 
scatter  in  the  yield  estimates  for  the  Russian  tests  is  thought  to  be  due  to  the 
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Figure  7.  Yield  estimates  for  the  period  1977-81  using  NTS  magnitude  versus  yield 
curves  corrected  for  differntlal  attenuation  between  NTS  and  Russian 
test  sites  (correction  applied  was  .35  magnitude  units  for  mb).  Yields 
based  on  Mg  were  obtained  using  "universal"  Mg  versus  yield  curves,  with 
no  correction  for  tectonic  release  effects. 


greater  perturbing  effects  of  tectonic  release  for  the  Kazakh  test  site  com¬ 
pared  to  NTS  or  the  Novaya  Zemlya  sites.  In  particular,  the  Kazakh  site  appears 
to  have  a  thrust  type  tectonic  release  equivalent  while  the  NTS  site  is  generally 
thought  to  produce  strike-slip  equivalent  radiation,  which  only  mildly  perturbs 
the  Ms  values  and  has  almost  no  effect  on  mb  values.  On  the  other  hand,  a 
thrust  mechanism  can  strongly  bias  I<is  values  and  have  nearly  as  large  an 
effect  on  mb  (increases  in  mb  would  be  expected).  As  observed  previously,  we 
intend  to  more  systematically  investigate  the  effects  of  tectonic  release  on 
these  yield  estimates  in  an  effort  to  provide  a  procedure  for  reduction  or  elimi¬ 
nation  of  its  effect  for  yield  estimation. 

V.  "Anomalous  Event"  Conditions  Leading  to  Discrimination  Failure  and  Inaccu¬ 
rate  Yield  Estimates 

It  is  rather  easy  to  demonstrate  that  there  are  several  complicating  effects 
that  can  cause  some  of  the  standard  earthquake-explosion  discrimination 
methods,  such  as  mb  vs  Ms,  to  fail  in  particular  circumstances.  Further,  biased 
yield  estimates  can  be  obtained  when  such  magnitude  data  is  used  to  estimate 
explosive  yield.  In  particular,  aside  from  the  obvious  failure  that  will  occur  for 
mb  vs  Ms  discrimination  when  earthquakes  are  grossly  misclassified  with 
respect  to  depth,  it  is  shown  that  particular  types  of  small  earthquakes  occur¬ 
ring  at  shallow'  depths,  in  the  range  from  5  to  25  km,  will  have  anomalously  low 
Rayleigh  wraves  in  the  usual  measurement  range  from  18  to  20  seconds.  This 
"Mt  null”  occurs  for  small  strike  slip  earthquakes,  and  is  found  to  occur  at  all 
azimuths  from  this  type  of  event.  Further,  as  illustrated  in  Figure  10  small 
thrust  earthquakes  with  dips  at  or  near  45D  produce  similar  "Ms  nulls”.  These 
low  Ms  values  will  lead  to  earthquakes  that  appear  very  ”explosion-like"  in  the 
m6  -M%  discrimination  plane. 


Yield  Determination  from  mb  8i  Ms  for  U.S.  Explosions 

(1968-1977) 


Figure  8.  Yield  estimates  using  both  and  Ms  for  NTS  underground  explosions 
in  the  period  1968-77,  prior  to  the  limited  test  ban  agreement. 
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A  second  complication  that  affects  spectral  discriminants  at  teleseismic 
distance  ranges,  and  affects  teleseismic  explosive  yield  estimation  as  well,  is  the 
highly  variable  effective  Q  observed  across  continents,  as  well  as  similar  varia¬ 
bility  across  trench  zones  and  oceanic  ridges.  For  example,  on  the  basis  of 
results  from  Butler  (1983),  which  are  summarized  in  Section  111,  large  variations 
in  the  observed  effective  Q  occur  within  the  continental  U.S.,  not  only  from  geo¬ 
logic  province  to  province,  but  also  within  these  provinces  on  a  rather  fine  spa¬ 
tial  scale.  While  part  of  these  observed  variations  are  undoubtedly  due  to  com¬ 
plex  velocity  structure,  which  produces  scattering  and  focusing  effects,  it  is 
certain  that  the  most  important  broad  scale  effect  is  due  to  upper  mantle 
absorption,  which  occurs  nearly  exclusively  within  the  upper  mantle  low  velo¬ 
city  zone.  As  is  well  known,  differences  in  the  anelastic  absorption  and 
"scattering  loss"  processes  between  test  sites  can  result  in  the  necessity  of 
recalibration  of  empirical  magnitude-yield  relations,  to  take  into  account  the 
differences  in  absorption  (the  "Q  bias")  between  test  sites,  as  was  discussed  in 
Section  IV.  This,  consequently,  requires  fairly  detailed  knowledge  of  the  upper 
mantle  characteristics  at  a  test  site,  or  reliance  on  other  (robust)  methods  of 
yield  estimation  which  are  unaffected  by  upper  mantle  absorption.  Thus,  for 
accurate  teleseismic  yield  estimation  using  "Q  corrected"  m 6,  it  would  be 
necessary  to  have  effective  Q  information  for  an  area  of  interest  that  was  at 
least  as  complete  as  that  now  available  for  the  continental  U.S. 

In  addition,  low  Q  zones  with  associated  high  absorption  can  cause  the  high 
frequency  spectra  of  explosions  and  earthquake  P  waves  to  be  reduced  below 
station  noise  levels  near  low  Q  regions,  so  that  the  spectral  differences  upon 
which  variable  frequency  spectral  discrimination  is  based  cn  be  obscured  by 
noise.  Such  effects  on  spectral  magnitudes  have  been  observed  and  correlated 
with  the  existence  of  anomalously  low  upper  mantle  Q  zones  in  the  vicinity  of 
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particular  events.  (In  this  case,  explosions  will  appear  within  the  earthquake 
population  in  the  spectral  magnitude  "parameter  space.") 

Theoretically  it  has  been  shown  (Stevens,  1980)  that  zones  of  high  initial 
prestress,  at  or  near  the  failure  zone  for  an  earthquake,  can  give  rise  to  peaked 
earthquake  spectra  and  a  significant  enrichment  of  the  high  frequency  part  of 
the  spectrum  from  the  event.  This  is  potentially  a  problem  for  both  mb  vs  Ms 
and  spectral  discrimination,  since  the  high  frequency  enrichment  could  lead  to 
enhanced  mb  values  at  a  fixed  Ms  level.  However,  it  appears  in  reality  that  the 
stress  levels  in  the  earth  at  shallow  depths  are  not  generally  high  enough  to 
cause  serious  difficulties  for  discrimination,  and  that  only  for  quite  deep  earth¬ 
quakes  (100  km  and  deeper)  is  the  effect  likely  to  be  very  strong  and  only  then, 
apparently,  for  some  rather  small  percentage  of  earthquakes.  In  this  regard, 
however,  Figures  11  and  12  illustrate  results  of  variable  frequency  magnitude 
(VFM)  discrimination  where  this  effect  appears  to  occur. 

Specifically,  these  figures  show  VFM  discrimination  results  at  Red  Lake. 
Ontario  and  the  LASA  array  in  Montana  for  well  recorded  events  from  a  large 
Eurasian  data  set.  At  RKON  three  deep  earthquakes  appear  in  the  explosion 
population,  plus  one  earthquake  that  is  thought  to  be  a  shallow  event. 
Nevertheless,  a  large  number  of  deep  and  shallow  earthquakes  form  a  distinctly 
separate  population  distribution,  which  is  well  separated  from  the  explosions. 
The  LASA  results,  with  fewer  events  well  recorded,  show  that  there  is  again  a 
good  representation  of  the  explosion  and  earthquake  populations  with,  how¬ 
ever,  one  deep  earthquake  appearing  in  the  explosion  population.  In  both  fig¬ 
ures  the  earthquake  populations  are  made  up  of  earthquakes  at  all  depths, 
from  very  shallow  to  very  deep.  Thus  it  appears  that  some  earthquakes  can 
violate  this  spectral  discrimination  criteria  at  a  single  station.  However,  it  is 
also  true  that  only  very  rarely  does  a  particular  earthquake  fail  the  VFM 


discrimination  criteria  over  an  entire  set  of  stations  and  in  the  one  instance 
where  this  occurred,  the  event  in  question  was  a  very  deep  earthquake  that 
appeared  in  explosion  populations  at  many  stations.  (The  event  was  in  fact 
event  No.  147  shown  in  these  figures.)  Thus,  there  appears  to  be  occasional 
deep  events  having  spectral  properties  that  are  explosion-like  and  that  appear 
so  at  nearly  all  teleseismic  stations.  Nevertheless,  almost  all  deep  and  shallow 
earthquakes  appear  distinctly  different  spectrally  from  explosions.  Clearly,  if 
these  "anomalous"  earthquakes  are  only  deep  earthquakes,  then  depth  esti¬ 
mates  for  all  events  that  appear  "explosion-like”  in  the  VFM  plane  will  identify 
these  "anomalous"  events  as  deep  earthquakes,  with  any  shallow  event  in  this 
population  being  a  probable  explosion. 

Should  there  be  more  common  occurrences  of  shallow  earthquakes 
appearing  in  the  explosion  population,  such  as  the  one  event  considered  to  be 
an  earthquake  in  the  explosion  population  region  in  the  VFM  plane  at  RKON. 
then  more  elaborate  (joint)  discrimination  procedures  would  be  required  to 
sort  out  the  probable  explosions,  from  events  in  the  explosion  population 
region  in  the  VFM  plane.  However,  based  on  the  Eurasian  events  studied  so  far. 
no  shallow  earthquake  that  was  "well"  recorded  (with  signal  to  noise  power 
larger  than  2)  at  several  stations  has  been  found  to  appear  explosion-like  at  a 
majority  of  the  stations.  Thus,  occasionally,  at  one  or  two  stations,  a  shallow 
earthquake  will  appear  explosion-like,  but  on  a  multi-station  identification 
basis,  these  events  have  been  properly  identified.  (This  means,  however,  that 
single  station  event  identification  by  VFM  spectral  discrimination  can  result  in 
misclassification  of  some  shallow  earthquakes  as  possible  explosions  unless 
other  single  station  discriminants  are  also  used  which  can  help  provide  identifi¬ 
cation.) 

Finally,  a  complicating  effect  of  considerable  significance  is  tectonic 


release  associated  with  explosions.  As  is  well  known,  numerous  Shagan  River 
test  site  explosions  show  highly  anomalous  surface  wave  radiation  involving 
large  Love  waves  and  "reversed"  Rayleigh  waves  (i.e.  reversed  in  polarity  com¬ 
pared  to  explosion  events  in  the  same  region  that  do  not  produce  large  Love 
waves).  The  effect  can  be  explained  by  tectonic  release  equivalent  to  a  thrust 
type  tectonic  shear  stress  field  at  the  explosion  site.  Hence,  it  is  clear  that  low 
frequency  Rayleigh  waves,  and  Ms  measurements,  can  be  very  strongly  affected. 
In  fact  the  tectonic  component  can  be  much  larger  than  the  explosion  com¬ 
ponent,  so  that  yield  estimation  from  Ms  can  be  (very)  strongly  perturbed.  In 
addition,  as  shown  by  Murphy,  Archambeau  and  Shaw  (1983),  the  direct  P  wave 
from  the  Rulison  explosion  in  Colorado  was  strongly  perturbed  by  induced  tec¬ 
tonic  release,  such  that  the  observed  mb  value  was  reduced  by  about  .3  magni¬ 
tude  units.  Therefore  both  Ms  and  mb  can  be  strongly  perturbed  by  tectonic 
release  effects  making  yield  estimation  using  these  magnitude  measures  more 
difficult,  (Indeed,  if  the  thrust  oriented  shear  stress  applies  to  the  surface 
waves  from  the  Shagan  River  events  then,  based  on  the  Rulison  results,  this  tec¬ 
tonic  mechanism  would  also  strongly  perturb  the  values  for  these  events. 
The  perturbation  would,  however,  be  such  as  to  increase  them  by  as  much  as 
several  tenths  of  a  magnitude  unit,  which  is  opposite  to  the  decrease  in  mb 
observed  for  Rulison.) 

Solutions  to  most  of  these  problems  are  not  too  difficult  to  find,  once  the 
mechanisms  responsible  for  the  anomalies  are  identified.  Specifically,  mul¬ 
tivariant  spectral  magnitude  discrimination  using  both  body  and  surface  wave 
magnitudes  measured  at  a  number  of  frequencies  can  eliminate  or  minimize  the 
discrimination  problems  described  earlier.  Here  the  appropriate  approach 
would  be  to  rely  upon  body  wave  spectral  magnitudes,  such  as  mb(f  ,)  vs 
mt,(f  g),  for  discrimination  in  the  regional  distance  range,  particularly  for  small 
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events,  and  to  use  mb(f  i)  vs  Ms(f  2)  (including  the  standard  mb  vs  Ms)  at 
teleseismic  distances,  for  the  larger  events  (eg.,  for  events  with  mb  >  4.0  to 
4.5).  In  this  case  the  "Mt  null"  for  small  events  is  avoided,  since  body  waves  are 
used  for  spectral  discrimination  of  small  events,  and  at  regional  distance 
ranges  the  upper  mantle  low  Q  absorption  will  not  affect  this  crustal  body  wave 
discrimination  procedure.  For  all  larger  events,  which  will  involve  failure  over  a 
fairly  wide  depth  interval,  the  "Ms  null”  will  not  appear  and  vs  Ms  will  apply 
as  a  well  verified  discrimination  procedure  that  is  not  so  strongly  perturbed,  by 
upper  mantle  Q  variations  and  other  propagational  effects,  as  to  cause  event 
identification  failure.  Here,  the  standard  vs  Ms  approach  would  be  supple¬ 
mented  by  using  mb  vs  Ms  at  a  range  of  frequencies  (e  g.,  1  Hz  and  .05  Hz,  as  is 
the  practice  now,  but  also  at  .5  Hz  and  .03  Hz,  1.5  Hz  and  .05  Hz,  etc.)  This  latter 
multivariant  mb  vs  Jlfs  approach  can  be  used  to  increase  the  confidence  level  of 
teleseismic  event  identification,  and  also  to  supplement  the  body  wave  spectral 
discrimination  at  regional  distances. 

For  accurate  yield  estimation  of  identified  explosions  in  all  situations,  it  is 
necessary  to  avoid  Q  biased  measurements,  or  to  correct  for  the  ”Q  bias",  and 
to  explicitly  remove  contamination  effects  due  to  tectonic  release,  the  latter 
when  it  occurs  being  evidenced  by  SH  wave  production  at  rather  high  levels  at 
all  frequencies.  This  can  be  accomplished  by  routinely  measuring  Ms(f )  (r  e  , 
Love  wave  magnitudes  at  several  frequencies)  as  well  as  M§(f)  (the  normal  Ray¬ 
leigh  wave  magnitudes  at  various  frequencies)  and  mb(f  ).  When  the  Love  wave 
magnitudes  for  identified  (or  suspected)  explosions  are  low,  then  M§(f ),  with 
appropriate  distance  corrections  applied,  can  be  used  to  estimate  yield  in  the 
usual  manner.  Since  based  yield  estimation  avoids  the  strong  upper  mantle 
attenuation  effects,  it  is  to  be  preferred  for  yield  estimation.  However,  mb 
based  yield  estimates  can  simultaneously  be  obtained  using  a  standard  empiri- 


cal  relationship,  and  any  "Q  bias"  for  the  explosion  site  region  can  be  estimated 
by  comparing  the  yield  estimated  from  m b  to  that  obtained  from  Ms.  (In  this 
way  a  calibration  for  mb  vs  yield  for  the  site  could  be  obtained.)  For  explosions 
showing  moderately  large  Mk  at  the  lower  frequencies  used  for  yield  estimation, 
such  as  .05  Hz,  then  the  observed  M§  and  Ms  data  must  be  jointly  interpreted 
(e.g.  within  the  theoretical  framework  of  stress  relaxation  due  to  shatter  zone 
induced  tectonic  release)  to  provide  estimates  of  the  physical  parameters  for 
the  tectonic  release  source,  in  particular  the  effective  shatter  zone  radius  and 
the  prestress  magnitude  and  orientation.  These  parameters  can  then  be  used, 
again  within  the  framework  of  the  theory,  to  infer  the  magnitude  of  the  explo¬ 
sive  component  for  both  M§  and  mb  measurements.  These  corrected  M§  and  mb 
values  can  then  be  used  to  infer  the  proper  yield,  again  with  M§  used  as  the  pri¬ 
mary  means  of  yield  estimation,  so  that  "Q  bias"  effects  are  minimized.  For  very 
large  observed  Mj;  values  associated  with  an  identified  explosion  (e  g.,  Mk  >  Ms 
for  /  ^  .05  Hz),  then  it  is  possible  that  earthquake  triggering  could  have  taken 
place  and,  in  this  rare  instance,  an  in-depth  study  of  such  a  potentially  complex 
event  would  be  necessary  to  obtain  an  accurate  explosive  yield  estimate. 
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I.  Mathematical  FonniatiaM 


In  the  following  sections  we  will  wish  to  make  use  of  a  number  of  rather 
standard  results  in  the  mathematical  theory  of  elasticity.  In  addition,  we  will 
use  some  new  results,  as  well  as  specialized  notation,  which  may  not  be  immedi¬ 
ately  familiar  to  the  reader.  In  order  to  present  the  more  physical  arguments 
and  discussion  that  will  follow,  without  the  clutter  involved  with  the  parentheti¬ 
cal  digressions  required  to  define  notation  and  develop  certain  mathematical 
results,  we  give  a  brief  summary  of  the  basic  mathematical  notation  and  results 
that  will  be  needed  later. 

Basic  Equations  of  Motion  and  Boundary  Conditions 

In  order  to  develop  the  representations  to  be  used  in  a  systematic  and  logi¬ 
cal  fashion  we  begin  with  the  basic  equations  of  motion  in  an  elastic/an  elastic 
continuum,  which  simply  are  an  expression  of  the  conservation  of  momentum. 
These  equations  can  be  put  into  the  form  (Archambeau  and  Minster,  197B): 

i«rur  =  Pt « ■•“•?=  1  *2.3.4  (1) 

where  the  summation  convention  for  repeated  indices  applies  and  where  ua  are 
components  of  the  displacement  field  and  fm  the  external  force  field  density, 
with  both  of  these  fields  written  as  four  vectors  for  the  convenience  of  express¬ 
ing  the  equations  of  motion  in  compact  form.  In  particular,  these  fields  are 
defined  to  have  the  component  form: 

(uT)  *  (u„  Ug,  u3.  0) 


(/.)  =  (/i./f*/s.  0) 


where  the  fourth  component  (the  "time-like"  component)  is  set  to  zero.  In 
addition,  the  independent  variable  vector  x  is  the  four  vector: 

(*7)  ■  (x,  .a*  *s.  O 


with  time  t  the  fourth  component  in  this  notation.  Further,  in  this  notation 
Is  the  "elastic  operator",  and  has  the  form: 


L** E  8x m  Ox, 


where 


= 


Cafyt  s  Glfltl i  (ot«  7<  s  (ii^i  t,  l) 
Gl4k4  ~  C««4 k  ~  P&ik'<  <*  *  =1,2,3 
Ctfyt  ~  0:  all  other  a,  (3,  y,  6 


Here  the  greek  indices  are  always  defined  to  run  over  the  integer  range  from  1 
to  4,  with  the  regular  latin  indices  i,  j,  etc.,  always  defined  to  run  over  the 
integer  range  from  1  to  3.  The  tensor  is  the  ordinary  elastic  tensor,  which 
in  the  isotropic  case  is  just: 

Qjki  *  Ad^dj*  +  t  +  dfl-dj i) 


(This  material  tensor  can  incorporate  anelastic  effects  if  the  Lamd  constants  X 
and  fi  are  generalized  to  represent  operators,  involving  convolutions  and 
derivatives  with  respect  to  time.  See  Archambeau  and  Minster  (197B)  for 
details.)  Finally,  p  is  the  density  for  the  material. 

Thus  the  elastic  tensor,  when  written  in  the  expanded  "four  component" 
form,  can  be  made  to  include  inertial  effects  as  well  as  rheological  effects.  The 
symmetry  of  this  generalized  elastic  tensor  is  the  same  as  that  of  the  ordinary 
elastic  tensor,  in  particular: 

QfiU  =  Cfikl  s  Cyu  -  Qan 

and 

=  =  = 

This  symmetry  insures  that  angular  momentum  is  conserved,  as  well  as  linear 
momentum,  in  the  equations  of  motion. 

A  generalized  stress  tensor  t ^  can  be  defined  in  terms  of  the  generalized 
elastic  tensor,  where 

duy 

=  te  7 

In  terms  of  this  tensor,  the  equations  of  motion  become 

T«*p*p/«.  (2) 

where  the  comma  denotes  the  derivative  with  respect  to  the  index  /?,  that  is: 

9 

As  is  shown  by  Archambeau  and  Minster  (1B7B),  boundary  conditions  at  or 
across  any  type  of  boundary  within  the  medium  are: 

I  7V7*  I  =  0; re  9K  (3) 

where  the  double  bracket  notation  is  used  to  denote  the  change  (or  jump)  in 
the  quantity  enclosed  across  the  internal  or  external  boundaries  BV.  That  is. 
[  p  ]  *  -p**\  where  p  is  any  function  or  variable  and  the  superscripts  (l) 

Bind  (2)  are  used  to  denote  the  value  of  p  when  the  surface  is  approached,  in 
the  limit  sense,  from  opposite  sides.  Here  r)»  are  the  components  of  a  space- 
time  normal  to  a  hypersurface  in  (x* ,  t ).  and  this  four-vector  has  the  form: 

(Vf)  s  (»i.  nf.  nB,  -  E$*n,) 

where  rq,  I  *  1,  2,  3.  is  the  ordinary  spatial  normal  to  a  surface  of  discontinuity 
in  the  continuum  and 

Vf=  Vi  -v,  . 

where  Vi  are  components  of  the  velocity  of  the  boundary4  and  vt  are  the 

*ffince  U  lathe  velocity  of  a  surface,  it  is  defined  to  be  continuous  across  that  surface. 

Thus  I  U\  I  s  0,  for  all  f . 
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components  of  the  Telocity  of  the  material  particles.  If  the  boundary  moves 
with  the  particles  and  the  particle  velocity  is  continuous  across  the  boundary, 
which  is  ordinarily  the  case  for  solid-solid  boundaries,  then  the  relative  velocity 
components  Up  all  vanish,  and  the  boundary  condition  (3)  reduces  to  the  regu¬ 
lar  condition  of  continuity  of  traction  across  a  material  boundary.  Further, 
since  the  combination  Upni  appears  in  the  expression  for  ijp,  then  Upni  van¬ 
ishes  if  only  the  normal  components  of  U  and  t  are  equal  and  continuous  at  the 
boundary.  This  is  the  case  for  ordinary  liquid-solid  boundaries.  However,  in 
certain  cases  such  as  when  the  boundary  in  question  is  a  moving  phase  change 
boundary,  where  the  material  on  one  side  is  in  a  different  phase  state  from  that 
on  the  other,  then  the  boundary  movement  in  space  may  be  driven  at  a  slower 
or  faster  rate  than  the  particles  are  driven,  and  in  this  case  Upni  is  non-zero. 
In  particular,  in  the  case  of  spontaneous  failure  in  a  solid  under  stress,  the 
boundary  of  the  failure  zone  expands  with  time  and  separates  material  regions 
in  which  the  rheological  properties  are  quite  different.  Further,  the  failure 
cone  surface  separates  regions  in  which  the  material  behaves  nonlinearly  (the 
interior)  from  the  exterior  linear  region,  where  the  linearized  equations  of 
motion  (1)  apply.  The  boundary  of  a  failure  zone  usually  expands  at  a  rate  that 
is  very  different  than  the  local  velocity  of  the  material  particles  and,  as  shown 
by  Archambeau  and  Minster  (1976)  and  in  later  sections  of  this  discussion,  the 
rate  of  movement  is  much  larger  than  the  particle  velocity,  usually  being  of.the 
order  of  the  local  shear  velocity  in  the  material.  Therefore  the  boundary  condi¬ 
tion  (3),  which  expresses  conservation  of  momentum  at  a  boundary,  must  be 
used  in  its  most  general  form  when  we  deal  with  failure  processes,  such  as  those 
involved  in  earthquake  phenomena. 

If  we  use  the  definition  for  the  generalized  stress  and  "elastic"  tensors, 
then  the  boundary  condition  (3)  can  be  put  into  the  form: 

I  WJ*  |  =  |  (pvkvp  -Tu)nt  ]  =  0 

where  Tit  is  the  ordinary  stress  tensor,  that  is 

Tit  —  CucnuiU nt,n 


in  the  exterior  linear  region,  and  where* 

vp  *  V|  -  Vi  *  -  Up 


Whereas  (3)  expresses  the  conservation  of  momentum  across  a  boundary, 
the  conservation  of  mass  is  insured  by  the  condition  (Archambeau  and  Minster, 
1978) 

I  pvp*i  1*0  (4) 


which  holds  in  general,  whether  the  boundary  moves  with  the  particles  or  not. 
If  the  boundary  does,  however,  move  with  the  particles,  then  conservation  of 
mass  is  insured  by 

[  v,nj  |  =  0 

*It  should  be  emphasised  that  U  ia,  by  definition,  continuous  scrota  the  boundary  over 
which  it  is  defined,  ia  that  |  Ui  J  =  0,  for  all  l .  However,  since  component*  of  V  cen  be 
discontinuous,  Ui  —  V|  can  only  occur  if  V  happens  to  be  continuous  in  al]  ita  components 
•n  the  surface  over  which  U  is  defined. 


•e  that  (4)  degenerate!  to  the  condition  of  continuity  of  the  normal  component 
of  the  particle  Telocity.  If  the  boundary  ii  assumed  to  be  a  "welded”  solid-solid 
contact,  then  the  tangential  components  of  the  particle  Telocity  across  the 
boundary  are  also  continuous,  and  the  (stronger)  condition, 

f  |-0.*  *1.2.3. 

applies.  This  is  also  often  expressed  in  terms  of  the  displacement  field  as 

[  ]  *  0.  *  *  I,  8. 3  . 

The  eonserration  of  energy  at  an  expanding  or  moving  boundary  has  the 
form  (Archambeau  and  Minster,  1978): 

|  (pEuf-vtTu  ♦  qt)nt  |  «0  (5) 

Here  p  is  the  density,  E  the  total  energy,  ft  the  components  of  the  heat  flux 
and  Tm  the  stress  tensor.  In  most  elastodynamic  applications,  including  those 
involving  failure  processes,  the  heat  flux  term  may  be  neglected  since  the 
changes  due  to  heat  flow  are  very  small  on  the  time  scale  associated  with  elas- 
to dynamic  variations.  Thus,  for  a  "fast  process"  like  failure  in  solids,  the  heat 
flux  term  in  (5)  is  negligible.  Further,  the  total  energy,  pE,  is  the  sum  of  the 
internal  energy  ps.  the  kinetic  energy  l/EptfeV*  and  the  potential  energy  due 
to  conservative  body  forces,  so  that 

JT*s*  4  1/EVftV*  <f  p 

where  p  is  the  potential  for  the  body  forces,  defined  by  /  *  —  Vp.  In  the  elasto¬ 
dynamic  case,  p  corresponds  to  the  gravitation  potential  and  changes  in  this 
quantity  across  boundaries,  in  particular  across  failure  cone  boundaries,  are 
very  small.  Therefore 

|£]ai|e|4  1/2  I  w*v*] 

is  a  good  approximation  for  the  change  in  E  across  material  boundaries.  When 
heat  exchange  and  changes  in  the  body  force  potential  are  neglected,  then  (5) 
is  automatically  satisfied  when  the  momentum  and  mass  conservation  relations. 
(3)  and  (4),  are  satisfied,  if  there  is  no  change  or  "jump"  in  the  internal  energy 
across  the  boundary.  This  is  the  situation  with  normal  boundaries,  defined  by 
material  contrasts  in  a  layered  medium,  where  the  boundary  moves  with  the 
material  particles. 

It  there  is  energy  absorption  or  release  at  a  boundary,  such  as  is  the  case 
for  a  failure  process  where  energy  is  absorbed  by  the  material  in  the  transition 
to  the  more  disordered  (failed)  state,  then  (5)  is  not  automatically  satisfied  and 
the  conservation  of  energy  condition  at  the  boundary  constrains  the  rate  at 
which  the  failure  process  can  proceed.  In  particular,  the  energy  absorbed  per 
unit  mass  (  s  )  is  characteristic  of  the  material  and  the  particular  transition 
process  which  is  occurring,  and  so  is  equivalent  to  a  quantity  analogous  to  a 
latent  heat.  For  failure  processes  in  solids  (which  may  or  may  not  involve  first 
order  phase  transitions)  the  energy  required  for  the  process,  per  unit  mass 
that  undergoes  failure,  is  denoted  by  L.  Hence,  we  set:  f  e  )  *  L,  and  the  boun¬ 
dary  condition  in  (5)  can  bs  put  in  the  forme 

P{1KL  4  1/2  |  vkvk  |)  (£$  -*/*>)*,  7i*n,  -  qtnt  J 
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Here  p^  and  vS1*  refer  to  the  material  before  transition,  that  is  to  the  material 
outside  the  failure  zone.  Ve  define  the  rupture  rate  Up  to  be  the  rate  at  which 
the  failure  surface  expands  in  the  direction  normal  to  itself,  measured  relative 
to  the  particle  velocity  in  the  material  outside  the  failure  zone.  Thus 

Up  «  (Vi 


Now,  using  the  other  conservation  relations  given  in  (3)  and  (4)  to  eliminate 
terms  involving  |  v*  ],  it  is  shown  by  Archambeau  and  Minster  (1976)  that  this 
relation  is  equivalent  to: 

« _  hi£i  _  jMj. 

w  pL  ZptL 

where  p  *  pi*!,  and  b  Tunt  is  the  traction.  Here,  of  course,  it  is  assumed  that 
L  t  0.  If  the  beat  flux  is  neglected  then  this  gives: 


which  expresses  energy  conservation  on  a  moving  boundary  along  which  energy 
absorption  occurs  due  to  a  transition  process,  such  as  material  failure. 

The  momentum  and  mass  conservation  relations,  in  the  general  case,  may 
be  combined  into  one  relation  by  expanding  (3),  and  by  using  (4)  to  eliminate 
terms.  Specifically,  the  combined  mass-momentum  boundary  relation  is 
(Archambeau  and  Minster,  1978): 

(?) 


where  Up,  p  and  tk  are  as  defined  above. 

In  general,  any  boundary  along  which  energy  absorption  (or  release) 
occurs  will  move  at  a  rate  that  is  different  than  the  particle  velocities  on  either 
side  of  the  boundary,  and  in  this  case  the  general  forms  of  the  conservation 
laws  expressed  by  (3),  (4)  and  (6),  or  (6)  and  (7),  are  the  appropriate  boundary 
relations  to  be  satisfied.  Further,  such  moving  boundaries  make  the  dynamical 
system  non-linear,  since  the  movement  or  expansion  of  these  boundaries 
depends  on  the  elastodynamic  field,  via  the  (nonlinear)  dependence  of  Ur  on 
the  particle  velocities  and  tractions  in  (5)  and  (7).  Thus,  while  the  equations  of 
motion  in  the  medium  exterior  to  a  failure  zone  are  linear,  and  are  as  indicated 
in  (l),  nevertheless  the  system  of  equations  governing  the  elastodynamic  fields 
in  the  linear  zone  outside  the  failure  region  consist  of  the  linear  differential 
equations  plus  the  (nonlinear)  boundary  conditions  on  a  rapidly  moving,  or 
expanding,  boundary  surface.  The  non-linear  aspect  of  this  elastodynamic 
problem  is  displayed  explicitly  by  the  Green's  function  integral  equation  for  the 
displacement  field  u  in  the  exterior  (linear)  region  surrounding  an  expanding 
failure  zone,  since  it  incorporates,  in  a  single  equation,  the  equation  of  motion 
and  boundary  conditions  at  the  failure  surface,  plus  initial  conditions. 

Green's  Functions  and  Integral  Equations  for  Oastodynamic  Problems 

In  order  to  obtain  integral  equations  for  the  elastodynamic  field,  and  exact 
and  approximate  solutions  of  these  latter  equations,  we  need  to  consider 
Green's  functions,  C^(x,  t ;  Xo>  *o)>  corresponding  to  the  k  component  of  dis¬ 
placement  in  the  medium,  at  a  receiver  point  at  x  and  time  t ,  due  to  an  impulse 


force  in  the  m  coordinate  direction,  occurring  at  a  source  point  at  Xo  and  at 
time  f0.  The  Green’s  function  to  be  considered  must  satisfy  the  inhomogeneous 
elastodynamic  equation  of  motion  along  with  the  regular  boundary  conditions 
for  a  layered  earth  model.  More  specifically,  the  boundaries  of  the  medium  are 
taken  to  move  with  the  material,  so  that  there  is  not  relative  movement  between 
particles  and  the  boundary  surface,  but  there  may  be  slip  at  the  boundary  such 
as  occurs  at  fluid-solid  or  fluid-fluid  boundaries.  In  this  latter  case  the  tangen¬ 
tial  components  of  the  particle  velocity  need  not  be  continuous. 

The  Green's  function  of  interest  therefore  is  defined  by  the  following  sys¬ 
tem  of  equations: 


tt  r 


'  afr 


-  AT(x.  t ixq.  t0) 


with  Ai"  a  generalized  delta  function  given  by 

AT  *  i(*  -  Xo)  «(f  ~*o) 

Here  6*n  si  if  i  =  m  and  is  zero  otherwise,  while  6(x  — Xo)  and  6(f  - 1 0)  are 
Dirac  delta  functions.  The  appropriate  boundary  conditions  for  CP  ore: 

I  9?  I  -  0 :  x  t  OV  (Bb) 

and,  on  welded  solid-solid  boundaries 

|CTl=;xzdV  (Be) 

or,  on  liquid-solid  and  liquid-liquid  boundaries 

|  sO-.xeBK  (Bd) 

m 

Here  (Jt  denotes  the  tractions  on  the  boundaries  of  the  medium,  the  latter 
denoted  by  OV.  with  n*  the  components  of  the  normal  to  such  a  surface.  (The 
boundaries  OV  can  be  internal  as  well  as  external,  with  any  of  these  being 
defined  as  surfaces  across  which  the  material  properties  are  discontinuous.) 
The  explicit  form  for  the  Green’s  traction  is: 


9*  *  w* 


The  system  of  equations  in  (B)  can  be  solved  in  closed  form  for  media  com¬ 
posed  of  layers  with  uniform  properties,  such  as  are  usually  used  to  model  the 
earth.  In  particular,  for  a  spherical  earth  model,  consisting  of  concentric 
spherical  shells,  the  Green's  function  satisfying  the  system  of  equations  (B)  in 
such  a  model  is  obtained  by  first  applying  the  Fourier  time  transform  to  the 
system  of  equations  (8),  and  then  expressing  the  Green's  function  as  an  eigen¬ 
function  expansion  with  coefficients  which  are  determined  by  the  boundary 
conditions.  The  resulting  system  of  algebraic  equations  can  be  solved  using  the 
Haskell-Thompson  matrix  method  (e.g.  Ben-Menabem  and  Singh,  1972)  and  this 
gives,  for  the  transform  of  the  Green's  function: 

«£)ff(xb:  «£) 

tt**  »■-(■%)« — 


.  f/(*  «£)  if(*o :»£) 


(9) 


where 

¥/(*  »£)  *  A»(r;  u£)  Ps»(d.  *)•* 

4  £**(r;  »&)  BW(*.  p)*^ 

f/(®  u£)  *  Ftm(r;  uD  Csn(d,  p)-s* 

Here  i>p  end  pj  denote  components  of  the  spheroidal  end  toroidal  eigenfunc¬ 
tions  for  the  layered  sphere,  with  and  the  associated  eigenfrequencies. 
A  bar  over  the  eigenfunction  denotes  complex  conjugation.  The  functions  Ptm, 
Btn  and  are  the  vector  spherical  harmonics  (t.g.  Horse  and  Feshbach.  1953) 
end  the  radial  functions  Dim,  E**  end  F^  can  be  expressed  in  closed  form  in 
terms  of  spherical  Bessel  functions,  with  coefficients  expressed  in  terms  of  pro¬ 
ducts  of  Haskell-Thompson  matrices,  (e.g.,  Ben-Henahem  and  Singh,  1972). 

An  entirely  similar  eigenfunction  expansion  for  the  Green’s  function  in  a 
layered  half  space  can  also  be  obtained  (Archambeau  and  Stevens,  unpublished 
manuscript,  1980).  Therefore  the  Green’s  functions  for  quite  detailed  models  of 
the  earth  can  be  expressed  explicitly  in  terms  of  elementaiy  functions. 

As  is  well  known,  the  system  of  equations  defining  the  Green's  function, 
that  is,  equations  (B).  may  be  used  together  with  the  elastodynamic  equations 
(1)  and  appropriate  boundary  conditions,  such  as  given  in  (8)  and  (7),  to  obtain 
an  integral  equation  for  the  unknown  elastodynamic  displacement  field  in  an 
elastic  (or  anelastic)  medium.  Such  an  integral  equation  incorporates  the  equa¬ 
tions  of  motion  along  with  the  spatial  boundary  conditions  and  "initial  values". 
This  integral  equation  is  usually  called  an  "elastodynamic  representation" 
theorem  in  the  present  context.  In  the  case  of  an  elastic/anelastic  solid  with 
an  expanding  failure  sone,  the  usual  development  of  the  Green's  function 
integral  equation  is  not  straightforward.  Archambeau  and  Minster  (1978)  give 
the  derivation  for  this  case.  The  important  differences  between  this  derivation 
and  the  more  standard  results  used  in  elasticity  theory  and  seismology  (e .g. 
DeHoop,  1957)  are  that  the  equilibrium  stress  field,  for  a  medium  which  is  in  a 
state  of  initial  stress,  will  change  as  an  internal  failure  sone  boundary  grows. 
This  gives  rise  to  an  "initial  value"  integral  term  that  isn’t  present  in  the  usual 
formulations  in  elastodynamics.  In  addition,  the  boundary  conditions  on  the 
growing  failure  surface  are  of  a  more  general  form  than  those  on  a  normal, 
"fixed",  boundary,  as  was  indicated  earlier. 

A  form  of  the  Green’s  function  integral  equation  which  accounts  for  these 
affects  is  shown  by  Archambeau  and  Minster  (1978)  to  be: 

4*um(m)  *  /  dt0  j p/«(xo)  (*•  *o>  *S*c  *o)  ty(*c)  *ao 


0Gf(x.*o)  .. 


(10) 


Here  the  summation  convention  applies  as  always,  and  a  four-vector  notation  is 
used,  where  all  greek  indices  run  over  the  range  from  1  to  4.  In  this  case  we 
express  ordinary  spatial  vectors  as  four  vectors  with  the  fourth  ("time-like") 
component  set  to  sero,  such  as  for  an  ordinary  displacement  vector  where: 

(«**)  *  (it  I.  U|.  v*.  0) 

Similarly,  all  ordinary  spatial  tensors  are  "adjoined"  with  zeros  for  the  "time¬ 
like"  (index  4)  components.  Thus: 

IGT\  if  a  or  p  *  1,2,3 
0;ifsorps4 

where  the  latin  indices,  such  as  1  and  m,  run  only  over  the  range  from  1  to  3. 

The  quantities  ijp  and  Jf,  on  the  other  hand,  are  four  vectors  with  non-zero 
time-like  components,  and  are  defined  by: 

Jf  - 

(ijp)  *  (fij,  nt.  »g.  — 

where* 

-  ,  * 

T«#  *  o)  ' 

,  ac^(x,x o) 

ffmf  *  £ap)rf(>b) 

are  the  generalized  stress  tensors  associated  with  the  displacement  and 
Green's  displacement  ffm  ,  respectively.  Both  and  ffmp  have  non-zero  "time¬ 
like"  components,  as  well  as  "space-like"  components.  Here  the  coordinate 
components  x«  define  a  four-vector  x  such  that 

(*s)  *  (*t.  *S.  **.«):«*  1 ,2,3,4 

and  Cmfr 4  Is  the  generalized  elastic  tensor,  defined  earlier  in  connection  with 
aquation  (1),  with  the  density  as  well  as  the  ordinary  elastic  constants  of  the 
material  as  components.  Thus  this  tensor  includes  inertial  as  well  as  elastic 
response  coefficients.  The  space-time  normal  with  components  rip  have  values 
equal  to  the  ordinary  components  of  a  the  spatial  (outer)  normal  to  the  moving 
surface,  plus  a  fourth,  "time-like",  component  with  a  value  equal  to  the  negative 
of  the  normal  component  of  the  boundary  velocity  at  this  surface. 

To  obtain  the  (linearized)  integral  relation  for  the  elastodynamic  field  u. 
the  equations  of  motion  in  the  linearized  form  of  (1)  are  used  and,  in  addition, 
the  generalized  normal  has  been  linearized  by  using  -  as  its  fourth 

*Wh#n  necessary  for  notation*!  clarity,  w*  vJU  write  the  source  coordinate*  a*  x  ® ,  xf ,  , 

*4  in  component  form.  Thaw  are  to  he  interpreted  a*  the  components  of  So-  We  shall  alto 
occasionally  writ*  So  as  Xr  when  necessary  for  consistency  Ja  notation. 


V  >  V  J* 
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component  rather  than  the  (proper)  term  -  Ufni,  with  Uf  ~  Vi  -  wj  where  vt  is 
a  component  of  the  particle  velocity  and  U\  a  component  of  the  (absolute)  velo¬ 
city  of  the  boundary.  In  the  case  of  interest,  wherein  we  are  concerned  with 
rupture  zone  growth  rates,  |U|  »  |  v|,  we  can  neglect  the  particle  velocity  rela¬ 
tive  to  the  rupture  velocity,  ao  Ufn*  *  C/jfij.  In  this  case  then,  the  components 
of  the  generalized  normal  rjf  can  be  expressed  in  a  form  that  is  explicitly 
Independent  of  the  elastodynamic  field  variables.  The  quantity  Ui  fit  is  usually 
called  the  "rupture  velocity",  and  is  defined  as  Ur  k  t/jnj .  It  is  simply  the  rate 
at  which  the  rupture  surface  expands  in  the  direction  normal  to  itself  and 
defines  the  time  history  of  expansion  or  growth  of  a  failure  zone. 

The  volume  integrals  in  (10)  are  defined  over  the  (linear)  region  txterior  to 
any  (nonlinear)  failure  zone,  and  this  exterior  volume  is  simply  denoted  by  V. 
The  surface  integral  is  to  be  taken  over  the  surface  0V  of  the  failure  zone  and 
all  other  boundary  surfaces  of  the  medium.  Because  the  failure  surface 
expands  with  time,  at  the  rate  Ur.  then  both  the  volume  V  and  the  surface  0V 
are  functions  of  the  source  time  variable  tg. 

• .  The  first  integral  in  (10)  involves  the  external  body  forces  acting  on  the 
medium.  In  the  elastodynamic  case  the  only  body  force  of  any  significance  is 
gravity,  and  only  the  changes  in  the  gravity  field  will  affect  the  dynamic  dis¬ 
placement  field.  These  changes  are  very  small  relative  to  tectonic  stress  field 
changes  for  an  earthquake  source  and  can  be  neglected  relative  to  other  terms 
in  (10).  (See,  for  example,  Archambeau  and  Scales,  19B3). 

The  final  integral  term  involves  the  derivative  of  the  equilibrium  displace¬ 
ment  field  u*  with  respect  to  the  source  time  variable  to*  When  a  failure  boun¬ 
dary  forms  and  expands  with  time  in  a  stressed  medium,  as  is  the  case  with  an 
earthquake  source,  then  the  equilibrium  state  for  the  medium  must  also 
change  and  this  integral  accounts  for  the  dynamical  effects  produced  by  such 
changes.  It  plays  the  role  of  the  initial  value  integral  of  the  classical  Green’s 
function  integral  representation  (e,p.  Morse  and  Feshbach,  1953).  but  has  a 
more  general  form  in  this  case.  The  surface  integral  over  the  moving  boundary 
OV  accounts  for  scattering  and  absorption  of  energy  at  the  moving  boundary 
and  it  also  has  a  more  general  form  than  arises  in  the  classical  case  due  to  the 
relative  movement  of  the  boundary  with  respect  to  the  material  particles. 

The  expanded  form  of  the  integrand  /Jfty  in  the  surface  integral  can  be 
shown,  from  the  definitions  given,  to  be: 


ft*,  58  GT  ^ 


Out  ) 


9*  ~P 


where  iq,  Tj*  and  fi*  are  the  conventional  spatial  displacement,  stress  and  sur¬ 
face  normal,  and  similarly  Q p*  and  ffu  are  the  Green’s  displacement  and  stress. 
This  quantity  is  identical  to  the  classical  result  for  a  stationary  boundary, 
except  for  the  addition  of  the  terms  involving  the  boundary  velocity  function  U. 
These  latter  terms  account  for  the  actual  transport  of  material  from  one  side  of 
the  boundary  to  the  other.  In  the  case  of  failure,  the  terms  account  for  the 
envelopment  of  new  material  by  the  rapidly  expanding  failure  zone. 

In  order  to  make  direct  use  in  (10)  of  the  eigenfunction  expansion  for  the 
Green’s  function  given  in  (9),  it  is  necessary  to  use  the  Fourier  transformed 
(spectral)  version  of  (10).  Hence,  with  the  Fourier  transform  operator  with 
respect  to  the  (observer  or  receiver)  time  variable  t  having  the  form: 


V/AVV.Va'.nIsV.V/.'.V.V 


SW!? 


//{*)•*«  dtmfiu) 


with  inverse 

W/(»)l  */?(») s“*  d/  *  /(f) 

where  o  =  2rr/ ,  then  application  to  a  causal  Green’s  function  (see,  for  example, 
Archambeau  and  Minster,  197B),  G£(k  Xo)  ■  Gf (xk,  x§;  f  -  f0),  gives: 

/i  Jsf  (** .  *o :  t  - 1 0)J  *  s'4*-0  S£  (xk ,  x\ :  «) 

Here  S£  denotes  the  Fourier  transform  of  G{  with  respect  to  t.  Since  the 
transform  is  with  respect  to  the  time  variable  i,  then  this  operator  commutes 
with  the  integrals  and  derivatives  taken  with  respect  to  x0  and  f0,  so  that  the 
operation  can  be  taken  inside  the  integrals  appearing  in  (10).  It  is  then  not  dif¬ 
ficult  to  show  that: 

4irffm(x*.  o)  *  /  s-4**  dt o  ^  fi/i  2 £*  dax0 

♦tqp2y,£fc)Jti*da0 

As  is  clear  from  the  expression  of  each  integral  term  in  (11),  the  terms 
correspond  to  Fourier  transforms  of  integrals  over  a  time  dependent  volume  or 
surface  integral.  Here  the  transform  turns  out  to  be  with  respect  to  the  source 
time  variable  t«,  rather  than  with  respect  to  the  receiver  time  t.  If  the  depen¬ 
dence  of  the  volume  V  and  surface  OV  on  to  can  be  neglected,  and  some  fixed 
limits  of  integration  used,  then  the  transform  operation  can  be  moved  inside 
the  spatial  integrals.  In  this  ease  the  result  (11)  simplifies  to  spatial  integrals 
Over  products  of  Fourier  transforms,  as  in  the  classical  result  with  stationary 
boundariea  However,  in  general,  the  result  (11)  must  be  used  essentially  as  it 
stands. 

One  formally  rigorous  approach  in  dealing  with  the  evaluation  of  (11)  in 
general  circumstances  is  to  (continuously)  map  the  surface  0K'  into  some  fixed 
surface  that  has  a  convenient  shape,  such  a  sphere,  and  to  thereby  produce 
spatial  Integrals  that  have  fixed  limits  with  respect  to  to-  In  this  case  all  the 
functions  appearing  in  the  integrand  are  also  mapped  into  a  new  (time  depen¬ 
dent)  "space"  and,  in  addition,  the  Jacobian  of  the  mapping  appears  in  the 
Integrand.  Then,  however,  the  Fourier  transform  can  be  moved  inside  the 
integrals  and  applied  directly  to  the  resulting  integrand.  In  this  case  a  variety 
of  approximations  are  possible  since  the  Jacobian  of  the  transformation  can  be 
expanded,  as  well  as  the  Green’s  functions  and  other  terms,  in  appropriate 
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function*  to  yield  quite  accurate  analytical  and  numerical  results.  This  iatter 
approach,  while  complicated,  has  been  pursued  by  Scales  and  Dilts  (1983)  and 
has  provided  initial  results  in  several  important  cases  Involving  ellipsoidal 
failure  cone  formation.  In  the  present  discussion  we  will  be  content  to  consider 
simpler  (geometric)  cases  that  can  be  approximated  without  recourse  to  the 
mapping  theory,  but  that  nevertheless  contain  the  essential  features  of  the 
theory  of  spontaneous  elastodynamic  sources  (earthquakes). 

In  view  of  the  length  of  the  expressions  for  the  integral  representation  of 
the  elastodynamic  field,  it  is  convenient  to  introduce  an  inner  product  notation, 
along  with  the  four-vector  (space-time)  representation,  in  order  to  condense 
the  results  and  to  highlight  their  physical  significance. 

Thus,  defining  inner  products  over  the  failure  zone  boundary,  0  V  and  the 
volume,  V,  exterior  to  that  boundary,  we  can  express  the  integral  representa¬ 
tion  of  the  transformed  dynamic  field  in  (11)  as: 

♦«&>  +  «£>  (12) 

with 

Sr)KJ 

»  iuF,,|(pa,Bu/.  cr>rJ 

Here  denotes  the  dynamical  displacements  associated  with  (gravitational) 
body  force  perturbations  due  to  the  failure  process.  As  noted  earlier,  this  term 
Is  very  small  compared  to  the  displacement  term  associated  with  the 
relaxation  of  tectonic  stress  due  to  the  spontaneous  creation  of  a  failure  zone. 
The  displacement  term  denotes  the  (secondary)  displacement  effects  due 
to  scattering  from  the  (growing)  failure  surface.  (The  four  vector  notation  is 
used  to  compress  the  expression  for  this  term.  However,  ujfi  is  identical  with 
vjf),  m  =  1,2,3  since  the  component  with  fi  *  4  is  zero.) 

The  various  inner  products  represent  spatial  integrals  over  the  volume 
outside  the  failure  zone  or  over  its  surface,  so  that  for  example: 

Sr),  •  ^0 

Finally,  Ftt  represents  the  Fourier  transform  operator  with  respect  to  the 
(source)  time  variable  1 0,  as  defined  earlier. 

If,  as  was  also  noted. earlier,  the  volume  V  and  surface  0V  are  independent 
of  time  variable  ic.  then  the  operator  Ft0  commutes  with  the  inner  product 
integrals  and,  in  this  special  case, 


«g>  ■  (fit..  §&,')„ -<r*».  8f>„. 

<tg> 

m* 

where  /j,  etc.,  ere  simple  Fourier  transforms  of  the  variables  in  question.  In 
this  ease,  when  the  volume  V  and  surface(s)  0  V  are  independent  of  time,  then 
it  is  also  true  that  the  normal  boundary  velocity  U-n  will  be  equal  to  the  normal 
component  of  the  particle  velocity  at  the  boundary.  That  is.  there  are  no  grow* 
ing  failure  sone  boundaries  In  this  ease.  Therefore  the  rupture  velocity 
Ug  =  (Ui  vanishes  in  the  integral  terms  above,  and  these  integrals 

immediately  reduce  to: 

*£>  «</>/*.  Sr  >r 

k^.Sr)^.  (w) 

where  Tu  and  (7"  are  the  ordinary  elastic  stresses  associated  with  the  displace¬ 
ment  14  and  Green’s  displacement  function  <**.  Since  the  medium  does  not,  in 
this  case,  have  an  internal  growing  Inclusion  .corresponding  to  a  failure  zone  (or 
other  phase  change),  then  the  equilibrium  field  will  not  change  with  the  time 
variable  f0  even  if  the  medium  is  initially  stressed.  In  this  case,  0iQu*  =  0.  and 
the  relaxation  displacement  field  vanishes.  Hence,  in  this  case  only  the  first 
two  integral  terms  in  (13)  are  nonzero  and  the  total  displacement  field  is  just 
the  sum  of  these  two  terms.  Written  out  explicitly  we  have,  in  this  "standard" 
oase: 

^  pfl  &r  «*8*c  + 

This  is.  in  fact,  the  classical  integral  representation  used  by  DeHoop  (1957)  and 
others  In  elastodynamic  and  seismological  applications.  As  will  be  seen,  this 
Integral  form  for  the  elastodynamic  problem  has  many  uses,  including  the 
representation  of  dynamical  dislocation  fields  and  equivalent  elastic  sources. 
In  these  cases,  and  in  many  others  involving  elastic  scattering  problems,  the 
first  integral  over  the. body  forces  ft  can  be  neglected.  On  the  other  hand,  if 
one  uses  equivalent  forces  to  represent  a  source,  then  this  leading  integral 
term  gives  the  displacement  field. 

II.  Basic  Theory  and  Models  of  Earthquakes 

It  is  generally  agreed  that  the  seismic  radiation  produced  by  an  earth¬ 
quake  is  a  consequence  of  the  relaxation  of  tectonic  stress  in  the  medium  sur¬ 
rounding  a  zone  of  material  failure.  Thus,  the  origin  of  the  energy  of  an  earth¬ 
quake  is  the  stored  strain  energy  within  a  large  volume  surrounding  the  (nar¬ 
row)  fracture  zone  and,  in  quantitative  terms,  the  total  energy  release  is  the 
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difference  between  the  equilibrium  etrain  energies  stored  in  the  entire  medium 
before  and  after  the  creation  of  the  failure  cone.  The  computation  of  this 
energy  change  is,  in  principal,  accomplished  by  simply  computing  the  strain 
energy  stored  in  the  prestressed  medium  with  an  inclusion  (the  failure  zone) 
and  subtracting  it  from  the  strain  energy  initially  stored  in  the  stressed 
medium  without  the  inclusion.  In  both  computations  the  prestress  state  is 
specified  by  requiring  the  stress  field,  with  or  without  tha  Inclusion,  to  assume 
a  prescribed  fixed  form  and  value  at  great  distances  from  the  origin,  where  the 
finite  failure  cone  is  centered.  (See  Landau  and  Ufshitz,  1989;  Archambeau, 
1968  and  1973  or  Eshelby,  1957  for  examples  of  this  computation  in  specific 
eases.) 

The  creation  of  a  failure,  when  viewed  in  this  way,  can  be  shown  to  cause  a 
net  reduction  of  the  strain  energy  stored  in  the  medium.  This  difference  in 
stored  (potential)  strain  energy  must  necessarily  be  partitioned  between 
seismic  radiation  and  the  energy  required  for  the  creation  of  the  failure  zone. 
Such  a  calculation,  involving  the  differences  between  two  equilibrium  states, 
does  not,  however,  tell  us  how  this  energy  is  partitioned  between  these  two 
forms,  nor  does  it  tell  us  how  the  energy  in  the  radiation  field  is  distributed 
spectrally.  That  is,  it  says  nothing  about  the  manner  in  which  the  system  goes, 
dynamically,  from  one  equilibrium  state  to  the  other. 

On  the  other  hand,  this  fundamental  energy  computation  forces  us  to 
recognize  that  the  energy  in  the  seismic  radiation  field  comes  from  the  release 
of  stored  energy  from  the  volume  surrounding  the  failure  zone.  Thus  an  earth¬ 
quake  is,  quite  literally,  a  volume  source  within  which,  or  at  the  boundaries  of 
which,  our  measurements  are  usually  made.  This  simple  observation,  in  itself, 
suggests  that  any  direct  dynamical  formulation  describing  the  seismic  radia¬ 
tion  from  an  earthquake  will  involve  a  volume  source  distributed  over  the 
medium  surrounding  a  growing  failure  zone  or  fracture.  This  follows,  of  course, 
because  the  energy  literally  comes  from  the  elastic  region  surrounding  the 
failure  zone  and  any  direct  dynamical  formulation  of  the  problem  designed  to 
predict  the  seismic  radiation  should  show  this  explicitly. 

lQnematical  Representations 

In  contrast  to  dynamical  representations  of  earthquakes,  which  must 
Involve  the  physics  of  the  failure  process  itself  as  well  as  stress  wave  radiation, 
kinematical  dislocation  sources  are  often  used  to  represent  an  earthquake. 
This  approach  has  considerable  intuitive  appeal  since  an  earthquake,  as  well  as 
changing  the  stress  field  in  the  medium,  also  produces  a  displacement  "disloca¬ 
tion",  or  offset  in  the  medium  across  the  failure  zone.  Since  a  failure  zone  is 
usually  very  thin,  often  only  a  few  millimeters  in  thickness,  then  the  zone  can 
easily  be  considered  to  be  essentially  a  plane  and  the  displacement  to  be 
discontinuous  across  this  plane.  Clearly,  one  may  conceive  of  an  earthquake  as 
being  equivalent  to  a  displacement  dislocation  taking  place  along  a  plane  (the 
"failure  plane")  with  an  "appropriate"  time  history.  In  this  case,  one  has  at 
ones  disposal  the  knowledge  that  if  a  displacement  is  specified  completely  in 
terms  of  its  space,  time  and  magnitude  history  on  such  a  plane,  and  if  the 
(usual)  boundary  condition  requiring  continuity  of  tractions  is  assumed,  then 
this  is  sufficient  to  prescribe  the  elastic  wave  field  everywhere  in  the  medium 
•urrounding  this  "dislocation  source".  This,  in  fact,  follows  quite  simply  from 
the  usual  Green's  function  solution  of  the  equations  of  motion  in  an  elastic 
solid  with  fixed  boundaries.  (See,  for  example,  Morse  and  Feshbach  (1953),  for 
basic  Green's  function  solutions;  or  Aki  and  Richards  (I960),  for  detailed  dislo¬ 
cation  solutions.)  In  particular,  the  displacement  field  (neglecting  body  forces) 
is  given,  in  general,  by*: 

*Th e  elimination  convention  is  to  b«  applied  to  repeated  indices  throughout,  so  for  ezenr 
hie,  Utfh  «  t/jTlj  4  C/gflg  +  Vans  Here  also,  we  use  t*  -  t  +  t,  where  t  is  the  re- 

•ehrer  time  and  t  >  0  ie  a  small  time  increment  which  merely  serves  to  avoid  any  ainaulari- 
tfeaofGT.  7 
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where  Cf  is  the  displacement  Green’s  function,  the  Green’s  function  trac¬ 
tion  on  the  surface(s)  E,  given  by 


pJW^o)  =  Qfuixc) 


0GT(z.fvo>fo) 


with  Cjft i  the  elastic  tensor,  and  where  4  end  u*  are  the  traction  and  displace¬ 
ment  components  of  the  dynamic  field.  On  a  dislocation  surface  E*.  where  the 
volume  enclosed  by  the  surface  is  allowed  to  vanish,  the  traction  component 
differences  across  such  a  (two-sided)  surface  vanish  by  the  assumption  of  con¬ 
tinuity  of  the  traction,  and  the  surface  integral  giving  the  displacement  field  in 
the  medium  reduces  to: 


4irum(x.t)  = 
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where  Au*  is  now  the  difference  in  displacement  on  the  tiro  sides  of  the  disloca¬ 
tion  surface  E*. 

It  is  no*  evident  that  if  the  surface  E*  is  specified  along  with  the  displace¬ 
ment  offset,  both  as  functions  of  time  and  spatial  position,  then  the  integral 
can  be  evaluated  and  the  field  is  given  (predicted)  at  all  receiver  positions  x 
and  for  all  time  t. 

It  is  just  as  clear  however,  that  all  that  is  accomplished  by  the  evaluation 
of  this  integral  representation  is  to  propagate  the  assumed  displacements  on 
the  surface  £4  to  other  points  in  the  medium.  Thus,  by  assuming  a  particular 
displacement  space-time  history  one  prescribes  a  kinematical  description  of 
the  failure  process.  The  problem  with  this  is  that  we  do  not  necessarily  know 
how  to  do  this,  and  even  if  we  are  fortunate  enough  to  be  able  to  guess  a 
proper,  or  nearly  proper  displacement  space-time  history  using  intuitive  con¬ 
cepts,  we  still  do  not  know  how  to  relate  such  an  assumed  displacement  func¬ 
tion  to  the  basic  physics  of  the  process.  Indeed,  this  inability  to  relate  the 
assumed  dislocation  displacement  to  the  fundamental  dynamical  variables  asso¬ 
ciated  with  an  earthquake  immediately  shows  up  in  the  computation  of  the 
energy  changes  due  to  a  dislocation  in  a  prestressed  medium  In  particular. 
Steketee  (1958)  demonstrated  that  the  strain  energy  change  in  the  medium 
due  to  creation  of  a  dislocation  is  independent  of  the  prestress  and  depends 
only  on  the  displacement  jump  and,  furthermore,  that  this  energy  change  is 
always  such  as  to  increase  the  strain  energy  of  the  already  stressed  medium. 
This  increase  is  precisely  the  opposite  of  what  is  required  for  an  allowed  spon¬ 
taneous  process.  Taken  together,  the  increase  in  strain  energy  and  the 
independence  from  the  initial  stresses  in  the  medium  simply  means  that  when 
displacements  are  arbitrarily  specified  on  a  boundary  or  dislocation  surface,  in 
a  manner  independent  of  the  initial  stresses,  then  work  is  done  on  the  system 
In  order  to  create  the  dislocation.  The  key  here  is  that  the  dislocation  displace¬ 
ment  is  imposed  arbitrarily  without  relation  to  the  initial  stress  state,  or 
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subsequent  time  changes  in  this  stress  state.  Further,  because  no  changes  in 
material  properties  occur  in  the  idealized  dislocation  procedure,  then  there  is 
no  reason  for  the  initial  stresses  to  adjust  and  the  dislocation  field  is  simply 
linearly  superposed  on  the  initial  field.  Consequently,  a  dislocation  representa¬ 
tion  for  a  spontaneous  energy  release  source,  such  as  an  earthquake,  is  funda¬ 
mentally  unrelated  to  the  basic  physical  processes  actually  occurring,  that  is. 
the  dislocation  displacement  is  not  related  to  the  dynamical  changes  of  the 
■tress  field  in  the  medium  within  the  representation  theory  provided  by  the 
dislocation  model.  (To  provide  such  a  relation  one  has  to  solve  the  dynamical 
problem  from  first  principles  and  then  deduce  from  the  solution  the  equivalent 
dislocation,  with  suitable  definitions  of  how  changes  in  the  energy  state  of  the 
system  are  to  be  interpreted.) 


Dynamical.  Relaxation  Source,  Representations 

While  a  dislocation  and  other  similar  ldnematical  representations  of  earth¬ 
quakes  can  certainly  be  useful  in  some  instances,  such  equivalents  are  not  very 
helpful  if  we  wish  to  describe  the  source  from  first  principles  in  terms  of  basic 
physical  parameters  and  variables.  We  shall  therefore  consider  the  description 
of  earthquakes  from  the  point  of  view  already  suggested,  namely  as  a  volume 
source  of  elastic  energy  which  is  released  as  a  consequence  of  failure.  Further, 
since  we  recognize  that  this  is  an  interactive  process  in  which  the  readjustment 
of  stress  in  the  medium  is  accomplished  dynamically  in  such  a  way  as  to 
strongly  influence  the  growth  of  the  failure  zone,  we  seek  a  formulation  of  the 
problem  that  incorporates  such  interactions. 

In  this  regard,  the  formulation  of  the  full  dynamical  problem  as  an  initial 
value  problem,  as  described  by  Arcbambeau  (1988)  or  Archambeau  and  Minster 
(1978)  for  example,  provides  a  representation  of  the  seismic  radiation  in  terms 
of  a  volume  integral  over  the  region  surrounding  a  growing  failure  zone.  This 
theoretical  representation,  which  is  now  termed  ‘relaxation  source  theory’, 
■hows  that  the  actual  (volume)  source  function  for  an  earthquake  corresponds 
to  changes  in  the  equilibrium  displacement  field  in  the  medium  surrounding  a 
failure  zone,  where  these  changes  arise  from  the  growth  of  a  new  boundary  (the 
failure  zone  boundary)  in  the  medium.  In  particular,  the  Green’s  function 
integral  representation  for  the  dynamical  description  of  both  the  failure  zone 
growth  and  the  associated  seismic  radiation  due  to  relaxation  of  the  equili¬ 
brium  stress  field  is,  from  Archambeau  and  Minster  (1978): 


4fnim(x.t) 


f  ** 


acp  eqr 


riido.  q 


(1) 

with  boundary  conditions  on  the  failure  zone  boundary,  denoted  by  8V  and  with 
normal  n,  given  by 

pV»\vt\ »-i«*i.  (conservation  of  mass  and  momentum)  (2) 


Vi  «  [f*  f*J/(2p*Z,)  ,  (conservation  of  energy) 


(3) 


where  Um  *  Uni  is  the  normal  component  of  the  velocity  of  rupture  boundary 
growth  and  L  la  the  energy  per  gram  required  to  cause  the  material  to  undergo 
failure.  The  latter  two  equations,  (2)  and  (3),  are  statements  of  conservation  of 
momentum  and  energy,  respectively,  at  the  moving  failure  sone  boundary  and 
the  bracket  notation  is  used  to  denote  the  changes  in  particle  velocity  com* 
ponents,  v*.  and  traction  components,  f*.  across  the  failure  zone  boundary. 

In  the  integral  representation  in  (1)  for  the  dynamic  displacement  field  u„ 
at  any  point  within  the  elastic  region  surrounding  the  failure  zone, 
fig*  si  Gg*(z.t;sQ,fo)  represents  the  Green’s  function  (impulse  response)  for  this 
elastic  medium,  with  x  and  t  denoting  the  receiver  coordinates  and  time,  while 
So  end  to  denote  the  independent  set  of  source  coordinates  and  time.  The  first 
integral  term,  involving  the  surface  integral  over  the  failure  zone  boundary, 
accounts  for  the  interaction  between  the  elastic  wave  field  and  the  growing 
failure  zone  boundary.  Ibis  interaction,  or  wavefield  "scattering",  is  actually 
coupled  to  the  boundary  conditions  of  (2)  and  (3),  since  the  rupture  velocity, 
Ug  =  £$ni,  is  explicitly  present  in  the  integrand  and  also  determines  (implicitly) 
the  limits  of  the  surface  integral,  since  the  surface  OV  changes  with  time  in  a 
manner  determined  by  the  function  Um- 

The  second  integral  term  over  the  entire  elastic  region.  V.  exterior  to  the 
failure  zone  involves  the  time  rate  of  change  of  the  equilibrium  displacement 
field  u*(xo.£o)<  measured  relative  to  the  initial  displacement  field  associated 
with  the  prestress  state  of  the  medium.  This  integral  therefore  accounts  for 
the  primary  radiation  due  to  relaxation  of  stress  in  the  medium  surrounding 
the  growing  failure  zone,  and  as  the  failure  zone  continues  to  grow  then  the 
source  density  term  [duV9*ol  the  integrand  will  be  a  non-zero  and  changing 
function  of  fo.  the  source  time.  This  integral  term  therefore  corresponds  to  a 
(generalized)  initial  value  contribution  to  the  radiation  field  which  arises  from 
the  continuous  change  in  equilibrium  that  is  required  during  the  creation  of  a 
failure  sone  inclusion. 

Special  Relaxation  Theory  Models  and  "Instantaneous"  Failure  Sources 

To  relate  this  integral  term  to  more  familiar  initial  value  problems  it  is 
instructive  to  consider  the  limiting  case  in  which  a  spherical  failure  zone  is 
created,  in  an  initially  stressed  medium,  at  a  rate  equal  to,  or  greater  than,  the 
intrinsic  compressional  velocity  in  the  medium.  This  situation,  although  some¬ 
what  of  an  idealization,  approximates  the  creation  of  the  roughly  spherical 
shatter  zone  by  an  underground  explosion.  In  this  case  the  shock  wave  pro¬ 
duced  by  the  explosion  propagates  at  a  speed  near  the  compressional  velocity 
in  the  solid  and  produces  a  zone  of  failure  in  the  material  at,  essentially,  a 
"supersonic"  rate.  If  the  medium  is  initially  stressed,  then  no  changes  in  the 
prestress  state  outside  this  failure  zone  can  occur  until  after  the  complete  for¬ 
mation  of  this  zone,  since  any  relaxation  of  stress  proceeds  (radially  outward) 
at  a  rate  which  is  less  than,  or  just  equal  to,  the  rate  at  which  the  failure  zone 
formation  occurs.  In  this  case  the  failure  inclusion  in  the  medium  is,  in  effect, 
created  instantaneously,  insofar  as  relaxation  of  any  prestress  is  concerned, 
since  its  speed  of  formation  is  greater  than  the  speed  with  which  information 
concerning  its  existence  can  be  propagated.  Thus,  an  initially  stressed  medium, 
which  had  been  in  equilibrium  suddenly  finds  itself  in  a  state  which  is  out  of 
equilibrium,  at  say  the  time  f 0  «=  0,  due  to  the  presence  of  a  newly  created 
internal  inclusion  which  has  altered  effective  elastic  properties.  In  this  case 
the  Initial  value  for  the  displacement  field  is: 

u*(xo.f o)  *  |u{0)(*o)  -  u(,)(xo)]//(f0) 


where  a*0'  is  the  initial  displacement  aisociated  with  the  prestressed  state  of 
the  medium  without  the  failure  zone  inclusion,  while  is  the  equilibrium  dis¬ 
placement  field  in  the  medium  with  the  failure  zone  present.  Further,  H(t )  is  a 
step  function,  at  f  •  =  0,  arising  from  the  "instantaneous"  creation  of  the  failure 
zone,  so  that  u*  is  zero  before  the  failure  zone  is  created  at  t0  =  0.  while  it  is 
constant  in  time  and  equal  to  [n^  —  u^]  after  the  zone  is  created.  According 
to  the  initial  value  integral  term  in  (1)  then,  the  source  density  term  0u**/  Bt 0 
will  be  a  delta  function  at  the  time  i0  *  0.  Thus,  the  initial  value  integral  term 
in  equation  (1)  takes  the  form 

W’telSlb- 
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where  K(0)  denotes  the  volume  integration  limits  at  the  time  f0  =  0,  and  so 
corresponds  to  the  entire  medium  outside  the  failure  zone.  For  all  times  f0  >  0 
the  rupture  boundary  velocity  Ur  is  zero,  since  the  failure  zone  is  fixed  in  size 
after  its  "instantaneous"  formation.  Therefore  the  boundary  conditions  in 
equations  (2)  and  (3),  for  this  special  problem,  reduce  to  the  (standard)  condi¬ 
tions  that  the  tractions,  f*.  are  continuous;  that  is,  to  M  =  0.  Likewise  the 
terms  in  the  surface  integral  involving  Ur  =  t$n(  are  absent;  so  for  this  "instan¬ 
taneous"  failure  problem  the  general  representation  of  the  displacement  field 
specified  by  equations  (1)  through  (3)  reduce  to: 

4m4m(*D,f)  ~fdt9  J  [c3"  f*  -u*  daD 
o  »rro) 1  1 


with  the  condition 


•cross  the  boundary,  OK,  of  the  failure  zone.  Here,  6u*  =  [u(0)(xo)  ~  u(I>(xo)] 
and 

Oil* 

tk  =  ***  oi/ n‘ 

are  the  factions  on  the  surface  with  normal  vector  n.  Similarly, 
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tre  Green's  function  tractions  on  the  surface  with  normal  nu 

The  equation  (4),  which  provides  a  representation  of  the  radiation  field  for 
this  special  relaxation  source,  is  identical  to  the  classical  Green's  function 
integral  representation  for  an  elastic  medium  with,  an  initial  displacement  from 
final  equilibrium  equal  to  6u*.  and  this  problem  is  discussed  in  standard  works 
in  considerable  detail  (s.p.,  Morse  and  Feshbach,  1953).  The  interpretation  of 
the  physical  origins  of  the  two  integral  terms  making  up  the  representation  in 
(4)  is  that  the  surface  integral  describes  scattering  of  the  primary  elastic  radi¬ 
ation  by  the  failure  sone  boundary  (*.g„  reflections  and  refraction  at  the  boun¬ 
dary),  while  the  second  integral,  involving  the  initial  value  of  the  displacement 
field,  represents  the  primary  elastic  field  produced  by  the  release  of  stored 
elastic  energy  from  the  entire  medium  external-  to  the  failure  xone.  Clearly, 
scattering  effects  will  be  of  secondary  importance  for  purposes  of  predicting 
the  main  features  of  the  seismic  radiation  field  from  this  relaxation  source  and 
to  first  order  one  expects  the  field  to  be  predicted  by  the  initial  value  term  in 
(4).  Therefore,  to  first  order,  we  have  the  simple  solution 

(6) 


where  denotes  a  first  order  solution  obtained  by  neglecting  the  surface 
integration  term  in  (4).  To  account,  approximately,  for  the  scattering  effects 
we  can  use  the  first  order  solution  for  ti*  and .f*  in  the  surface  integral  term  in 
equation  (4),  to  obtain  a  second  order  approximate  solution.  That  is.  we  have  as 
an  approximate  solution  which  includes  scattering  effects: 

4 mi£>(z.f )  *  4rruJk,,(x.O  +  f  dt0  [g5T*jF,)  -  “k(1)^r)rfa0  (7) 

The  detailed  solution  of  (8)  for  a  spherical  failure  xone  can  be  found  in  Arch  am- 
beau  (1972),  for  example,  and  involves  first  specifying  a  prestress  state  for  the 
medium  and  then  solving  the  static  boundary  value  problem  for  an  inclusion,  in 
the  prestressed  medium.  This  requires  a  solution  of  the  equations  of  equili¬ 
brium  such  that  the  new  equilibrium  stress  assumes  the  value  of  the  prestress 
at  large  distances  from  the  failure  xone  while  satisfying  continuity  of  traction 
and  normal  displacement  across  the  failure  sone  boundary.  Here,  in  particular, 
the  material  inside  the  failure  sone  is  considered  to  have  rheological  properties 
that  are  quite  different  than  those  before  failure.  Solutions  for  static  inclusion 
problems,  for  various  boundary  shapes,  can  be  found  in  the  literature  (e.p.,  Lan¬ 
dau  and  lifshitx,  1959;  Neuber,  1946;  Esbelby,  1957)  or  can  be  obtained  from 
first  principals  in  either  analytical  or  numerical  form.  In  any  case,  using  the 
static  solution  so  obtained  to  specify  the  initial  value  du*  in  (8),  as  the  differ¬ 
ence  in  the  equilibrium  states  before  and  after  the  failure  xone  formation,  and 
then  expressing  the  Green's  function  as  an  eigenfunction  expansion  (i.e.,  as  an 
expansion  in  eigenfunctions  for  a  layered  spherical  or  flat-earth  model,  for 
example),  allows  the  integral  in  (8)  to  be  evaluated.  Once  is  obtained  from 
(8),  then  the  second  order  solution  uj?*  can  be  obtained  from  (7),  using  the 
expansion  for  the  Green's  function  and  the  same  integration  procedures  as  are 
used  for  (8). 

It  is  appropriate  to  emphasise  that  the  integral  representation  for  the 
radiation  field  given  by  equation  (4),  along  with,  the  iterative  form  of  solution 
given  by  (6)  and  (7),  can  be  used  to  describe  the  seismic  radiation  from  any 
inotantanoouM  failure  process  with  any  boundary  shape.  That  is,  while  a 
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■pberical  failure  tone  ia  one  that  la  of  practical  intereat  and  can  effectively  be 
produced  inatantaneoualy  by  an  explosive  ahock  wave,  aucb  a  failure  zone 
abape  la  only  one  of  an  infinite  set  of  geometric  possibilities  that  can  be 
described  by  (4).  Thus,  while  it  ia  difficult  to  envision  natural  situations 
(including  man  made)  where  failure  zones,  other  than  explosive  induced  spheri- 
-  cal  (or  near  spherical),  are  produced  at  supersonic  rates  (or  "instantaneously") 
nevertheless  in  the  abstract  we  can  use  the  "geometric  generality"  of  the  sim¬ 
ple  "instantaneous  failure"  representation  in  other  ways,  as  we  will  show  later, 
to  provide  an  understanding  and  interpretation  of  the  rather  complex  general 
representation  in  (l). 

Interpretation  of  the  General  Relaxation  Source  Representation 

The  value  of  the  "instantaneous  failure"  example  is  at  least  two  fold:  First 
it  shows  that  the  rather  complicated  representation  given  in  (l),  for  a  general 
spontaneous  failure  process,  reduces  in  a  natural  and  logical  way  to  what  is 
essentially  a  classical  result  with  which  we  are  familiar  and  for  which  there  are 
alternate  formulations  that  can  be  demonstrated  to  give  the  same  result.  In 
this  regard  Stevens  (1961),  has  shown  that  the  initial  value  formulation  given  by 
(4)  can  be  put  into  the  form  of  a  stress  pulse  equivalent,  wherein  the  initial 
value  volume  integral  can  be  re-expressed  as  a  surface  integral  over  the  failure 
zone  boundary,  and  where  the  effective  source  function  in  this  integral  is  now 
the  stress  difference  between  the  initial  and  final  equilibrium  states  at  the 
boundary.  This  form  is  identical  to  the  usual  integral  representation  employed 
in  crack  theory.  Further,  when  reformulated  in  this  way,  Stevens  was  able  to 
solve  for  the  seismic  field  produced  by  a  spherical  inclusion  exactly,  including 
all  scattering  effects.  Further,  his  solution  was  obtained  for  an  arbitrary  (spa¬ 
tially  variable)  prestress  condition,  so  that  the  solution  has  considerable  gen¬ 
erality.  Comparison  of  this  exact  solution  with  the  iterative  solution  procedure 
Indicated  by  the  equations  (6)  and  (7),  but  with.  (7)  repeated  in  successive 
approximations,  showed  that  the  iterative  method  is  convergent  and  that  the 
scattering  term  produced  only  a  small  correction:  to  the  relaxation  term,  as 
expected. 

In  addition  to  the  verification  of  the  general  representation  and  its  itera¬ 
tive  solution  In  a  special  case,  the  "instantaneous  failure"  example  can  be  used 
to  provide  insight  into  the  structure  and  meaning  of  the  general  relaxation 
source  representation.  That  is,  the  limiting  process  and  logic  used  to  generate 
the  special  case  equations  from  the  general  representation  of  (i)  through  (3) 
can  be  turned  around  to  enable  us  to  see  how  the  general  formulas  can  be  built 
from  this  special  case.  In  particular,  we  can  adopt  the  special  instantaneous 
failure  integral  representation  and  associated  solution  as  the  fundamental  or 
canonical  problem  to  be  used  to  build  more  general  solutions.  From  this  point 
of  view  the  more  general  case  of  spontaneous  failure  at  a  finite  rate  can  be 
thought  of  as  a  sequence  of  instantaneous  failure  processes  of  the  canonical 
type,  with  each  increment  of  failure  zone  growth  occurring  instantaneously 
after  some  interval  of  time  6f0  between  successive  instantaneous  events.  This 
corresponds  to  the  approximation  of  a  continuous  growth  rate  by  a  sequence  of 
■mall  steps.  If  we  add  up,  or  superpose,  the  fields,  produced  by  each  step-like 
change  in  the  failure  zone,  then  we  will  be  able  to  obtain  the  total  radiation 
field  from  the  entire  sequence.  Proceeding  with  this  approach,  using  (6)  to 
represent  the  radiation  induced  by  each  increment  of  boundary  growth  and 
considering  scattering  (and  other  interactions)  at  the  boundary  as  a  second 
order  effect,  we  get  to  first  order: 

*2^ p[*“.*CO] <,s*0 
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where  1,2,  •  •  •  N,  denote  the  times  at  which  incremental,  instantaneous, 

growth  of  the  failure  sone  occurs.  Here  V(tn)  denotes  the  elastic  volume  out- 
mida  the  failure  sone  at  the  time  t*.  Now  we  can  explicitly  introduce  the  time  of 
separation  between  increments  of  instantaneous  growth  of  the  failure  by 
observing  that  t m  ♦  6tc  and  rewrite  the  previous  superposition  as: 


d**e 


Now  we  can  pass  to  the  limit  in  which  dfe  i*  infinitesimal  and  where  the  boun¬ 
dary  growth  also  becomes  infinitesimal,  and  with  it  the  incremental  changes  in 
the  equilibrium  field  6uk*.  By  this  process  we  can  simulate  a  continuous,  and 
finite  rate,  failure  process  and  we  have,  using  the  definition  of  an  integral  as 
the  limit  of  the  summation: 

We  observe  now  that  this  limit  of  a  sequence  of  elementary  initial  value 
representations  of  the  seismic  radiation  from  instantaneous  failure  processes 
results  in  a  total  field  representation  that  is  precisely  the  same  as  the  relaxa¬ 
tion  term  in  (1).  which  accounts  for  stress  relaxation  effects  in  the  general 
case.  Thus  we  arrive  at  an  interpretation  of  the  most  important  term  in  the 
general  representation  in  terms  of  an  elementary  initial  value  source  which  has 
a  well  known  solution,  in  particular  the  general  source  is  just  a  "sum"  of  these 
elementary  initial  value  sources. 

The  close  relationship  between  the  very  simple  "instantaneous  failure" 
source  and  the  more  complicated  integral  representation  for  the  general  case, 
given  in  (1).  can  be  exploited  further  by  observing  that  the  surface  integral 
term  represents  scattering,  or  interaction  of  the- primary  field  with  this  boun¬ 
dary  discontinuity  in  both  cases;  albeit  that  the.  boundary  interaction  in  the 
general  case  is  more  complex  since  it  involves  reflections,  etc.  from  a  moving 
boundary  and  simultaneous  absorption  of  energy  at  this  boundary,  due  to  the 
energy  required  to  produce  a  failure  transition  in  the  material.  In  spite  of 
these  differences  however,  we  may  approach  the  solution  of  the  integral  equa¬ 
tion  in  (l)  in  precisely  the  same  fashion  as  was  done  in  the  case  of  the  simple 
"instantaneous  failure",  and  so  in  analogy  with  the  iterative  approach  leading 
to  equations  (6)  and  (7),  we  use,  in  the  general  case,  the  solution  procedure: 

4mdr>(«.l)  =  [«■  [«««  -£?-*>  ~iir~  V,\ 
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where  the  index  n  fc  2  is  the  order  of  the  iteration  carried  out  using  (9),  with 


(B)  the  first  order  "starting  approximation".  Here  the  exact  solution  is,  for¬ 
mally: 

«■»(*.<)  “  Jm 

but  in  practice  or  even  u^l\  is  a  sufficiently  good  approximation  for  the 
seismic  radiation  field,  just  as  was  the  case  for  the  instantaneous  failure. 
Hence  the  primary  contribution  to  the  radiation  field  arises  from  the  general¬ 
ised  initial  value  term  given  by  while  corrections  accounting  tor  scattering 
and  the  energy  adsorption  required  to  drive  the  failure  process  are  obtained 
using  the  iterative  approximation  prescribed  by  (9).  The  energy  absorption 
manifests  itself  explicitly  through  the  dependence  of  the  integrals  in  (B)  and  (9) 
on  the  time  varying  limits  of  integration,  which  are  determined  by  the  rupture 
rate  function  Ifa,  and  through  the  rupture  rate  dependent  terms  of  the  form 
pvt  Ug,  with  v*  =  0ut  /  8t  corresponding  to  the  particle  velocity,  appearing  in 
the  surface  integral  in  (9).  The  interpretation  placed  on  these  effects,  which 
were  not  present  in  the  instantaneous  failure  representation,  is  that  there  is 
strain  energy  absorption  in  addition  to  scattering  from  a  moving  failure  zone 
boundary  in  this  more  complete  description.  Further,  we  do  not  prescribe  the 
boundary  movement,  but  instead  the  rupture  rate  is  determined  by  the  dynami¬ 
cal  boundary  conditions  at  the  failure  surface,  given  by  equations  (2)  and  (3). 
In  this  interpretation  then,  energy  absorption  arises  because  the  dynamical 
radiation  field  itself  "feeds"  the  failure  process  and  supplies  the  necessary 
energy  to  create  a  larger  failure  zone.  Thus,  the  strain  energy  is  redistributed 
dynamically  by  wave  propagation  to  the  failure  sone  boundary  and,  in  one  or 
more  locations,  reaches  levels  high  enough  to  supply  the  energy  required  to 
continue  the  failure  growth  process.  When  the  failure  transition  occurs  (albeit 
in  a  very  narrow  zone,  but  nevertheless  in  some  finite  volume  of  the  material) 
then  the  local  strain  energy  is  absorbed  in  the  process,  essentially  being  used 
to  convert  the  material  to  its  "failed  state”. 

In  spite  of  the  added  complications  associated  with  the  interaction  of  the 
radiation  field  with  the  (nonlinear)  failure  process  in  the  spontaneous  failure 
representation,  the  approach  to  the  solution  of  the  integral  representation  for 
the  radiation  field  is  the  same  as  for  the  instantaneous  failure.  That  is,  the  first 
order  solution  is  obtained  by  first  solving  the  static  boundary  value  prob¬ 
lem  for  an  inclusion  of  appropriate  geometrical  shape  and  rheological  proper¬ 
ties  in  an  initially  prestressed  medium,  so  that  u**,  the  change  in  the  equili¬ 
brium  field,  is  obtained.  This  static  field  is  expressed  in  a  form  that  depends 
parametrically  on  the  time,  f0,  because  of  the  dependence  of  the  inclusion 
geometry  (i.s.,  rupture  zone  dimensions)  on  the  failure  velocity  function.  Thus, 
for  example,  the  equilibrium  field  changes  due  to  the  creation  of  an  ellipsoidal 
failure  zone  in  a  prestressed  medium  would  depend  parametrically  on  the  ellip¬ 
soidal  axes  (a,  b  and  c)  of  the  Inclusion,  and  these  in  turn  would  depend  on  the 
time  integral  of  the  rupture  velocity  in  the  axial  directions,  so  that  the  change 
in  the  equilibrium  field  would  be  expressible  in  terms  of  a  time  variable  t0.  The 
derivative  of  the  equilibrium  field  with  respect  to  this  time  parameter  yields  the 
required  source  term  in  (B),  and  with  the  Green’s  function  expressed  in  terms 
of  an  eigenfunction  expansion,  then  evaluation  of  the  integral  for  is  possi¬ 
ble.  Once  wJt*)  has  been  obtained,  in  the  analytic  form  of  an  eigenfunction 
expansion,  then  the  iterative  procedure  expressed  by  (9),  while  algebraically 
cumbersome,  is  straightforward. 
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bnplicatl«DS  from  the  Dynamical  Boundary  Conditions  for  Failure  Growth 

In  order  to  solve  this  source  problem  using  the  approach  described,  or  any 
other  method  for  the  solution  of  (1),  it  is  clear  that  a  determination  or  specifi¬ 
cation  of  the  rupture  velocity  function  Cfc(xo.to)  is  necessary.  As  was  indicated 
earlier,  the  boundary  equations  (2)  and  (3)  serve  as  dynamical  conditions  for 
failure  and  therefore  actually  determine  the  rupture  velocity  function.  In  most 
modeling  of  earthquakes  however,  the  boundary  conditions  have  not  been  used 
explicitly  to  determine  this  function,  but  instead  a  rupture  velocity  function 
has  been  specified  (or  assumed)  having  a  form  that  was  based  on  a  combination 
of  simple  observation  and  elementary  analysis  of  the  failure  boundary  condi¬ 
tions.  (See,  for  example,  Archambeau,  1988  and  Minster,  1973). 

Arcbambeau  and  Minster  (1978)  show  formally,  however,  that  if  a  complete 
loss  of  shear  strength  and  no  essential  change  in  compressibility  characterizes 
the  material  after  failure,  which  is  a  rather  plausible  assumption,  then  the 
boundary  relations  in  (2)  and  (3)  give  a  relation  of  the  form: 


where  vk  are  particle  velocity  components  and  C,  denotes  a  signal  velocity 
characterizing  stress  wave  propagation  in  the  medium  surrounding  the  failure 
surface.  If  the  rigidity  inside  the  failure  zone  is  essentially  zero  and  if  the 
failure  process  is  considered  to  be  driven  by  a  shear  wave,  so  that  C,  at  Vp/p, 
then  this  relation  gives: 


Ur  t*  'JyJ  p 


with  fj,  and  p  denoting  the  shear  modulus  and  density  in  the  medium  outside  the 
failure  zone.4  Therefore,  under  the  reasonable  assumption  of  loss  of  shear 
strength  in  the  material  involved  in  the  failure  process,  it  is  found  that  the  rup¬ 
ture  rate  will  be  close  to  the  local  shear  velocity  within  the  medium.  This  result 
appears  quite  consistent  with  observations  which  suggest  that  the  highest  rup¬ 
ture  rates  are.  in  fact,  close  to  the  shear  velocity. 

Another  result  that  provides  insight  into  the  failure  dynamics,  as 
expressed  by  the  conservation  equations  at  the  failure  boundary,  can  be 
obtained  from  the  energy  relation  (3)  which  states  that 

<12> 


Now  we  note  that  the  previous  result  states  that  Ur  t*  vt  »  Vp/p  when  the  pro¬ 
cess  of  failure,  characterized  by  a  particular  energy  of  transition  L0,  results  in 
a  complete  loss  of  shear  strength.  Thus  for  this  process 


e  change  in  compressibility  is  allowed  during  the  failure  transition,  then  thia  result  is 
modified  and  it  is  possible  for  the  rupture  rate  to  become  "trans-sonic”,  that  is,  to  have  a 
value  between  the  Shear  and  compress) ona]  velocities  in  the  medium,  ft  is  felt  that  transi¬ 
tions  of  thl*  sort  are  unlikely,  or  at  least  rare  and  confined  to  failure  transitions  at  great 
depth  in  the  earth. 


which  follows  from  the  general  energy  equation,  where  we  have  used  (Ac|  to 
represent  the  quantity  V|f*f*|  which,  in  itself,  is  a  measure  of  the  magnitude 
of  stress  drop  across  the  failure  boundary.  Because  of  the  lack  of  shear 
strength  in  the  failure  sons  for  this  process,  then  |Ac|  represents  a  "total 
stress  drop"  (that  is,  the  maximum  possible  shear  stress  change)  at  the  boun¬ 
dary. 

Now  for  some  other  process  of  failure  that  may  operate  in  some  thermo¬ 
dynamic  environments  in  the  earth,  a  complete  loss  of  shear  strength  may  not 
occur,  or  the  material  may  be  characterized  by  plastic  or  viscous  properties 
allowing  the  retention  of  shear  tractions.  In  this  case  the  energy  equation  still 
applies  in  the  form  given  by  (12).  with  Ur  and  the  traction  jump  in  general  dif¬ 
ferent  from  those  in  (13).  But  we  can  use  (13)  as  a  means  of  normalizing  the 
general  case  relation  given  by  (12),  so  as  to  obtain  a  scaling  law  relating  "par¬ 
tial  shear  stress  drops"  to  the  condition  of  total  stress  drop,  characterized  by 
| Ao|.  That  is,  if  we  divide  equation  (12)  on  both  sides  by  the  result  in  (13)  we 
get: 


14, 
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Thus  we  find  that  the  rupture  rate  scales  as  the  ratio  of  the  partial  to  total 
(shear)  stress  drops,  and  is  always  directly  proportional  to  the  local  shear  velo¬ 
city  of  the  material  surrounding  the  failure  zone.  We  also  observe  that  Ur  will 
depend  on  the  square  root  of  the  transition  energy  ratios  when  the  process  of 
failure,  resulting  in  a  partial  stress  drop,  is  fundamentally  different  than  the 
failure  process  resulting  in  complete  loss  of  shear  strength.  However,  in  this 
regard,  we  might  reasonably  make  the  assumption  that  the  same  fundamental 
process  is  normally  responsible  for  all  occurrences  of  failure  in  the  earth,  and 
that  the  differences  in  stress  drop  are  due  to  relatively  subtle  changes  in 
mineralogy,  water  content  and  local  temperature,  which  result  in  different 
rheological  characteristics  among  similar  earth  materials  subject  to  the  same 
process  of  failure.  In  this  case  we  would  have  L  *  L*.  and  transition  energy 
differences  would  not  mediate  the  normalized  rupture  rate  equation,  given  by 
(14).  Under  this  assumption  then,  we  get  the  simple  scaling  relation 

|  Ac.  | 
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where  we  have  introduced  JAc,  |  to  represent  the  magnitude  of  the  "partial 
stress  drop"  quantity  |f*f*  f • 

These  results  allow  us  to  make  several  generalizations  concerning  earth¬ 
quake  dynamics  and  to  also  simplify  the  analytical  modeling  of  these  sources  of 
seismic  waves.  In  particular,  total  loss  of  shear  strength  in  the  material  after 
failure  implies  total  stress  drops  and  rupture  rates  at  or  near  the  local  shear 
velocity.  The  converse  of  this  result,  namely  that  rupture  rates  near  or  at  the 
■bear  velocity  necessarily  require  total  loss  of  shear  strength  (and  total  shear 
■tress  •  drops)  has  not  been  shown,  but  certainly  this  would  be  a  plausible 
Interpretation  of  high  rupture  rate  observations.  Further,  it  is  found  that  rup¬ 
ture  velocities  Ur  scale  linearly  with  the  magnitude  of  the  stress  drop,  so  that 
the  larger  the  partial  stress  drop,  the  targer  the  rupture  rate,  until  the  limits  of 
total  stress  drop  and  rupture  rates  equal  to  the  shear  velocity  of  the  medium 
are  reached. 


-  We  observe  further,  however,  that  the  stress  drop  and  traction  jumps 
appearing  in  all  the  relations  describing  conditions,  at  the  failure  boundary,  in 
particular,  equations  (14)  and  (IS),  are  dynamical  quantities  and  actually 
correspond  to  transient  tractions,  or  stresses,  at  the  points  on  the  failure 
boundary  surface  where  active  failure  is  taking  place.  These  stress  changes 
.  may  be  very  large  due  to  dynamical  stress  concentration  effects  in  the  vicinity 
'?  of  the  expanding  edge  of  the  failure  cone  and,  while  clearly  related  to  the 
'  ambient  background  tectonic  stress,  could  be  one  or  even  two  orders  of  magni¬ 
tude  larger  than  the  ambient  stress  level.  Therefore  to  relate  an  equation  like 
(15)  to  ambient  quasi-static  tectonic  stresses,  it  is  necessary  to  assume  that 
the  transient  dynamic  stress  concentrations  occurring  at  the  front  of  a  failure 
sone,  which  are  responsible  for  failure,  are  directly  related  to  the  (local) 
ambient  stress  level.  The  simplest  assumption  would  be  that  the  general  shape 
and  curvature  of  the  failure  surface  is  the  same  for  all  conditions  of  failure,  in 
which  case  the  dynamic  stress  levels  are  always  (for  all  failure  modes)  directly 
proportional  to  the  ambient  stress  levels.  With  such  an  assumption,  the  ratio  of 
the  partial  to  total  dynamic  stress  drops  is  equal  to  the  ratio  of  the  changes  in 
the  ambient  (quasi-static)  stress  levels  for  partial  and  total  stress  drop.  In  this 
case  a  relation  like  (15)  can  be  used  to  estimate  ambient  tectonic  stress  levels, 
since  it  can  be  rewritten  in  the  form 

Vr  j“j“  (16) 


where  |Aff|  is  the  quasi-static  change  in  the  ambient  stress  field  under  condi¬ 
tions  of  a  partial  stress  drop  and  |  o«  |  denotes  the  ambient  stress  level,  which 
arises  from  the  total  stress  drop  condition.  .  •  • 


Physical  Processes  and  Stress  Level  Requirements  for  Failure  In  the  Earth 

If  assumptions  required  to  arrive  at  the  result  given  in  (16)  are  appropri¬ 
ate,  in  particular  if  the  mechanism  of  failure  has  the  same  energy  requirements 
whether  the  stress  drop  is  total  or  partial  and  if  the  dynamical  stress  load  caus¬ 
ing  failure  is  strictly  proportional  to  the  quasi-static  ambient  stress  for  any 
failure  process,  then  the  relation  provides  the  means  of  estimating  the  ambient 
stress  from  simultaneous  observations  of  the  rupture  rate  and  the  seismically 
observed  stress  drop.  Such  estimates  cannot,  however,  be  literally  interpreted 
as  being  the  stress  levels  that  directly  cause  failure,  since  it  is  the  dynamic 
stress  concentrations  that  produce  failure  once  the  process  is  initialed  and, 
quite  certainly,  it  is  intense  quasi-static  stress  concentrations  that  initiate 
failure  in  the  first  place.  That  this  must  be  true  follows  from  an  elementary 
consideration  of  the  energy  requirements  for  reasonable  processes  of  failure  in 
the  earth. 


Specifically,  we  observe  that  direct  and  indirect  estimates  of  tbe  average 
stress  change  in  the  quasi-static  stress  field  due  to  earthquakes  are  in  the 
range  near  100  bars  (10*  dynes/cm*)  and  are,  at  most,  from  300  to  500  bars 
(eg.  Archambeau  et  oi„  1963).  ftirther,  the  highest  rupture  rates  observed  are 
near  the  shear  velocity,  so,  using  (16),  we  estimate  that  ambient  deviatoric 
stresses  are  at  most  of  the  order  of  500  bars,  and  usually  closer  to  100  bars  in 
regions  with  strong  earthquake  activity.  Now  if  stress  levels  of  such  low  magni¬ 
tude  actually  were  directly  responsible  for  failure,  then  from  (13).  we  would 
have  for  the  energy  per  unit  volume,  pLo.  required  for  failure,  the  estimate: 


.  i  itoi»  l  (s « 

,L,“  8  fZr*  8  3(3  x  10s)* 


4.6  *10*  ergs /cm9 
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But  this  implies  an  energy  for  failure  in  the  earth  that  is  extremely  low  in  com¬ 
parison  with  energy  requirements  for  even  trivial  ohanges  in  the  thermo¬ 
dynamic  state  of  the  material.  For  example,  the  energy  required  to  raise  the 
temperature  of  one  (ram  of  water  by  one  degree  centigrade,  in.  say.  a  water 
■aturated  rock,  is  ltr  ergs,  which  is  already  25  times  greater  than  the  failure 
energy  estimated  from  the  largest  ambient  stress  field  changes.  Further,  when 
we  consider  the  energy  required  for  a  material  phase  change  such  as  shear 
melting,  then  even  in  the  presence  of  water  which  would  reduce  the  melting 
point  considerably,  the  latent  heats  Involved  are  of  the  order  of  10*  ergs/cm* 
and  larger,  and  so  the  energy  required  for  failure  by  melting  would  be  of  the 
order  of,  at  least,  10*  ergs/c to*.  This  energy  is  nearly  four  orders  of  magnitude 
larger  than  the  failure  energy  estimate  based  on  ambient  stress  field  changes. 

Similarly,  if  we  consider  the  ambient  strain  energy.  W.,  stored  in  the 
medium  per  unit  volume,  and  take  the  change  in  this  quantity  after  failure  to 
be  an  estimate  of  the  energy  per  unit  volume  available  for  initiation  of  failure, 
then  if  the  average  ambient  stresses  determine  failure  energies,  we  would  have 
to  have 

Here,  as  previously,  pL0  is  the  energy  per  unit  volume  required  for  failure  and 

is  the  average  tectonic  stress.  Hence,  with  the  average  ambient  stress 
taken  to  be  500  bars  as  an  extreme,  we  have,  for  a  total  loss  of  shear  strength 
at  the  failure  boundary: 

pL0  **•  |  IT* |  at  5  x  10”  erge/cm* 

Thus,  with  this  quite  different  method  of  estimating  the  energy  for  failure,  we 
obtain  essentially  the  same,  very  low,  transition  energy  estimate  when  reason¬ 
able  estimates  of  the  largest  average  ambient  stresses  are  used. 

Clearly,  one  conclusion  that  must  be  made  is  that  the  failure  energy 
inferred,  using  reasonable  upper  limits  for  ambient'  stress  field  changes  result¬ 
ing  from  earthquakes,  is  so  small  that  it  cannot  be  a  valid  estimate  for  the 
energy  required  for  failure  in  the  earth.  Nevertheless,  it  is  even  more  certain 
that  tectonic  stresses  must  be  responsible  for  failure  in  the  earth  and  so  we  are 
forced  to  conclude  that  stress  concentration  effects,  both  transient,  during 
failure  growth,  and  quasi-static,  prior  to  failure  initiation,  are  the  relevant 
stress  levels  supplying,  the  very  local,  high  energy  levels  required  for  failure. 

Indeed,  in  the  dynamical  situation,  where  we  are  concerned  about  the 
energy  required  for  continued  failure  growth,  the  failure  cone  already  exists 
and  it  can  be  viewed  as  an  inclusion  in  the  medium  having  drastically  different 
elastic  properties  than  the  surrounding  rock.  Such  an  inclusion  will  produce 
very  large  stress  concentrations  in  the  surrounding  medium  in  the  near  vicinity 
of  that  part  of  its  surface  having  large  curvature,  for  example  near  the  edges  of 
an  ellipsoidaliy  shaped  failure  inclusion.  For  a  growing  inclusion  corresponding 
to  the  expanding  failure  sone,  the  stress  concentrations  are  dynamically 
changing  and  transient,  but  certainly  very  large  at  the  ends  of  a  long  thin 
failure  transition  sone  with  a  width  of  a  few  millimeters.  Indeed,  for  such  an 
Inclusion,  stress  concentration  factors  of  10  to  100  could  be  expected.  In  this 
case,  the  initial  tectonic  energy,  which  as  we  have  seen  has  a  low  value  per  unit 
volume,  is  focused  by  the  failure  sone  itself.  Here,  the  strain  energy  is  dynami¬ 
cally  transferred  through  the  propagation  of  stress  waves  from  other  parts  of 
the  medium  to  the  ends  of  the  failure  sone.  We  know  that  this  spatial 


redistribution  of  internal  energy  is  assured  because  the  system  is  totally 
governed  by  the  necessity  of  maintaining  dynamical  equilibrium  with  the  inclu¬ 
sion  zone  and  so  must  satisfy  boundary  conditions  at  the  failure  inclusion  sur¬ 
face  that  demand  stress  distributions  that  approach  those  for  the  static  limit, 
these  being  characterized  by  large  stress  levels  at  inclusion  edges.  Therefore, 
the  process  of  stress  relaxation,  acts  to  increase  the  stress  levels  and  strain 
energy  levels  near  the  high  curvature  edges  of  the  inclusion  and  to  decrease 
the  stress  and  strain  energy  elsewhere,  with  the  very  local  and  large  energy 
increases  at  the  inclusion  edges  being  accounted  for  by  part  of  the  energy 
reduction  in  the  rest  of  the  medium  and  with  the  "excess”  strain  energy  reduc¬ 
tion  being  radiated  to  the  far  field.  Now  the  energy  momentarily  stored  at  the 
edges  of  the  inclusion  must  be  large  enough  to  supply  the  energy  required  for 
failure  and  when  this  occurs  the  failure  zone  will  grow. 

The  dvnamical  condition  for  failure  growth  at  a  rate  Ur  is.  in  fact,  given  by 
equation  (12).  so  that  we  can  estimate  the  magnitude  of  the  transient  stress 
concentration  required  to  continue  the  failure  process  at  some  specified  rate. 
In  particular,  if  we  consider  failure  progressing  at  or  near  the  shear  velocity  in 
the  medium  with  total  loss  of  shear  strength  occurring,  then  the  magnitude  of 
the  dynamical  (shear)  stress  drop  |  Ao|  must  be  given  by 

|  Acs  1  *  =  pv,  V2Z,0  (17) 

where,  as  in  equation  (13),  we  have  used  |Ao|  to  represent  the  square  root  of 
'  the  traction  jump  across  the  failure  surface. 

Since  the  failure  process  requiring  teast  energy  will  be  initiated  first  and 
will  prevent  stress  concentrations  adequate  to  trigger  other  failure  processes 
requiring  higher  activation  energies,  we  must  consider  the  process  with  the 
lowest  energy  value,  L0,  in  order  to  estimate  |Acr| .  In  this  regard,  it  seems  that 
the  lowest  energy  process  that  has  been  proposed  as  a  mechanism  leading  to 
failure  is  a  temperature  increase  in  a  water  saturated  rock  medium  sufficient 
to  drive  the  local  effective  stress  to  zero.  Under  the  supposition  of  very  low 
tensile  strength  for  the  material,  then  only  a  small  tensional  stress  would  result 
in  failure.  This  phenomenological  mechanism,  proposed  by  Raleigh  and 
Evemden  (1962)  on  the  basis  of  quite  a  different  set  of  arguments  than  those 
considered  here,  would  require  the  presence  of  free  water  and  very  low  permea¬ 
bility  in  the  rock.  However,  both  conditions  are  likely  for  the  material  in  tec¬ 
tonically  active  regions.  (Arguments  for  free  water  and  high  pore  pressure 
based  on  geochemical  and  other  evidence  are  given,  for  example,  by  Fyfe  c t  a L, 
1978.) 

For  this  phenomenological  mechanism  to  operate  it  is  necessary  to  also 
propose  a  secondary  mechanism  that  would  rapidly  convert  the  transient  elas¬ 
tic  strain  energy  at  the  edge  of  the  failure  zone  to  heat  and  produce  the 
required  increase  in  water  temperature.  Such  a  mechanism  would  quite  cer¬ 
tainly  have  to  be  a  microscopic  process,  activated  by  the  high  dynamic  stress 
levels.  Numerous  mechanisms  are  actually  possible,  particularly  when  we 
recognize  that  the  material  is  extremely  heterogeneous  on  a  microscopic  level 
end,  therefore,  would  be  expected  to  contain  microscopic  stresses  from  10  to 
100  times  the  local  average  stress,  which  is  already  at  a  much  higher  level  than 
the  ambient  tectonic  stress  (by  another  factor  of  from  10  to  100)  due  to  the 
"stress  focusing"  at  the  edge  of  the  failure  zone.  Thus  we  can  expect  extremely 
high  stresses  at  the  edges  of  grains  and  pores  which  will  activate  highly  dissipa¬ 
tive  flow  and  cracking  in  the  material  on  the  microscopic  level  In  fact,  it  would 
be  likely,  or  at  least  possible,  that  the  material  could  locally  melt  if  water  wasn't 


*60 


-re¬ 


present,  but  in  the  presence  of  free  water  then  the  heat  produced  by  the  non¬ 
linear  and  highly  dissipative  processes  would  rapidly  ruse  the  fluid  tempera¬ 
ture  and  with  it  the  fluid  pressure.  During  this  process  the  high  microscopic 
•train  energy  levels,  and  the  average  local  strain  energy  level,  would  drop 
rapidly  as  non-linear  plastic  flow  relaxed  the  stress.  When  the  local  interstitial 
fluid  pressures  reached  the  level  of  the  maximum  principal  stress  in  the  rock 
matrix  or  exceeded  it,  then  the  already  locally  weakened  and  microscopically 
flowing  material  could  catastrophically  yield,  or  fail,  with  the  localized  zones  of 
flow  and  microfracture  connecting  throughout  an  element  of  material  volume. 

This  particular  process  of  failure  would  have  a  relatively  low  energy  of 
activation,  that  is  a  low  failure  transition  energy  Zo-  In  this  regard,  Raleigh  and 
Evernden  (1982)  estimate  that  the  pressure  increase  required  in  the  water 
would  be  of  the  order  of  the  average  ambient  tectonic  stress,  since  the  water 
would  be  initially  at  a  pressure  about  equal  to  the  least  principal  stress  in  the 
rock  matrix.  Therefore  a  pressure  increase  of  the  order  of  100  bars  would  be 
required.  This  would  occur  with  about  5aC  Increase  In  temperature.  Now  if  the 
transient  strain  energy  stored  in  an  element  of  volume  at  the  edge  of  the 
failure  zone  Is  essentially  all  converted  to  heat  by  the  microscopic  processes  of 
plastic  flow  and  microfracture  (with  associated  stress  relaxation  due  to  this 
deformation  and  flow)  then  the  amount  of  energy  required  to  raise  the  tem¬ 
perature  by  5aC  in  a  water  saturated  rock  matrix  is  about  equal  to  the  total 
energy  required  for  failure,  that  is,  equal  to  Zo  in  our  notation.  Thus,  for  this 
process 

Zo  at  (5aC)  (4  x  10*  ergs/ eel)  (.25  co l/gm/'C)  *  5  xiO7  ergs/gm 

where  .25  cal/gm/*C  is  a  representative  specific  heat  for  the  predominantly 
rock  medium.  Certainly  this  is  a  low  energy  when  compared  to  the  energy  per 
gram  required  for  melting,  which  is  of  the  order  of  10*  ergs/gm  or  larger. 

Now,  given  that  this  process  requires  an  energy.  Zo  of  about  5  x  107 
ergs/gm,  tbeu  we  can  use  the  relation  (17),  which  expresses  energy  conserva¬ 
tion  at  the  expanding  failure  surface,  to  obtain  an  estimate  of  the  magnitude  of 
the  stress  concentration  required  for  the  failure  process.  That  is,  we  require 

|  Ao|  =  pvt  V2Z0  ***  3(3  *  I®8)  ^10®  *  9  kb 

Thus,  when  the  ambient  tectonic  stress  is  of  the  order  of  100  bars,  then  stress 
concentration  factors  of  the  order  of  100  are  required  in  order  to  achieve  the 
required  transient  stress  levels  in  the  lower  kilobar  range.  This  can  only  be 
achieved  with  thin  failure  zones,  so  that  even  this  relatively  low  energy  process 
of  failure  demands  stress  concentrations  of  such  high  levels,  relative  to  the 
ambient  stress,  that  only  thin  failure  zones  are  predicted  to  be  possible.  Such 
thin  transition  zones  are,  of  course,  what  are  observed  for  earthquakes,  where 
new  fault  zones  are  commonly  observed  to  be  of  the  order  of  millimeters  in 
thickness.  Therefore,  in  this  regard  at  least,  these  inferences  are  consistent 
with  observation.  We  also  note  that  the  inferred  kilobar  stress  levels  are  of  the 
same  order  as  are  observed  In  rock  mechanics  experiments  for  failure  at  high 
dynamic  strain  rates. 

The  probable  situation  that  must  occur  prior  to  the  initiation  of  failure  is 
likely  to  be  similar  to  the  description  of  the  dynamical  growth  of  the  failure 
zone.  In  this  case,  however,  there  is  no  pre-existing  failure  surface  and  no 
(rapid)  dynamical  stress  readjustments  taking  place  to  concentrate  the  stress. 
However,  it  is  clear  that  some  process  of  stress  concentration  is  necessary  in 
order  to  initiate  failure  in  view  of  the  previous  energy  considerations.  In  this 
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case,  with  quasi-static  stresses  involved,  it  is  likely  that  much  lower  stress  lev¬ 
els  are  sufficient  to  activate  non-linear  creep  and  flow  processes  on  a  micros¬ 
copic  level,  particularly  when  water  is  present.  Therefore  pre-existing  hetero¬ 
geneities  within  the  medium,  on  a  macroscopic  scale,  could  serve  to  concen¬ 
trate  the  very  slowly  increasing  average  tectonic  stress,  eventually  to  a  level 
required  for  activation  of  creep  and  micro-scale  fracture  processes.  Since 
these  mechanisms  would  be  expected  to  be  activated  in  a  heterogeneous 
manner  spatially,  then  the  resulting  spatially  non-uniform  local  weakening  and 
flow  of  the  material  could  be  expected  to  increase  the  stress  levels  even  further 
at  those  points  where  the  material  was  stronger  and  more  resistant  to  such 
deformation.  At  such  points  stress  levels  could  reach  higher  levels  on  a  shorter 
time  scale  and  new  microscopic  processes,  requiring  higher  stress  levels,  could 
be  activated.  At  some  stage,  with  the  assumed  presence  of  stronger  zones 
allowing  the  eventual  concentration  of  quite  high  stress  levels,  one  would 
expect  that  a  new  increment  of  load  accumulation  within  the  highly  stressed 
strong  cones  would  exceed  the  threshold  for  activation  of  a  new  set  of  micros¬ 
copic  deformation  mechanisms.  This  could  trigger  rapid  and  intense  deforma¬ 
tion  and  weakening,  producing  a  small  very  weak  zone  which  would  serve  as  a 
stress  concentrating  inclusion  on  a  macroscopic  scale.  Once  this  occurred 
then  the  dynamical  process  described  earlier  could  serve  to  continue  the 
growth  of  such  a  zone,  sustained  by  the  elastic  energy  released  from  the  sur¬ 
rounding  medium  and  concentrated  at  the  inclusion  edges.  In  this  view  of  the 
initiation  and  growth  of  a  failure  zone  then,  it  might  be  expected  that  there 
would  commonly  be  a  "starting  time  period"  during  which  the  failure  growth  was 
relatively  slow.  This  time  period  of  initiation  could,  however,  be  quite  short  in 
duration,  with  dynamical  stress  concentrations  rapidly  driving  the  failure 
growth  rate  to  the  shear  velocity  level. 

As  a  final  observation,  it  is  appropriate  to  point  out  that  the  rupture  rate, 
Ur,  can  be  expected  to  be  highly  variable  spatially.  Thus,  for  example,  while  the 
shear  velocity  appearing  in  the  equations  specifying  Ur  can  be  nearly  constant 
spatially,  the  traction  or  stress  magnitude  jumps,  and  the  ratios  for  partial  and 
total  stress  drop  magnitudes,  can  be  expected  to  vary  in  a  manner  related  to 
the  heterogeneous  spatial  variations  of  the  ambient  stress  field.  In  addition, 
the  partial  stress  drop  magnitude  might  be  expected  to  vary  in  a  manner 
dependent  on  relatively  subtle  material  properties  (e.g„  water  content  and/or 
mineralogy).  Therefore  these  results  suggest  that  while  the  dynamics  of  the 
spontaneous  failure  processes  resulting  in  earthquakes  obey  some  (superfi¬ 
cially)  simple  basic  relations,  these  same  relations,  when  combined  with  our 
knowledge  of  the  heterogeneous  makeup  of  the  medium,  predict  that  the 
failure  process  can  be  very  complex,  with  highly  variable  rupture  rates  and 
stress  changes  not  only  possible,  but  quite  likely  to  occur.  Therefore  we 
expect,  and  in  fact  observe,  that  large  earthquakes  involving  large  spatial  rup¬ 
ture  cones  will  be  very  complex  and  produce  seismic  radiation  fields  that 
reflect  highly  variable  stress  drops  and  rupture  rates  along  the  failure  zone. 
Indeed  many  large  events  look,  seismically,  like  multiple  events  with  "bursts"  of 
seismic  radiation  observed  from  different  areas  of  the  failure  zone  as  it 
expands  with  variable  speed.  Further,  it  is  likely  that  small  events  would  have 
the  same  character,  but  with  such  extreme  variations  taking  place  on  a  smaller 
spatial  scale  and  thus  on  a  shorter  time  scale,  so  that  these  effects  would  be 
seen  in  the  seismic  radiation  at  higher  frequencies,  and  likely  to  be  observed 
only  in  the  near  or  regional  distance  ranges  from  the  event. 


—a--  . 

t* 


Relaxation  Theory  Models  for  Earthquakes 

The  results  relating  the  rupture  rate  to  the  intrinsic  shear  velocity  in  the 
medium  and  the  stress  drops  occurring  during  failure  have  significant  implica¬ 
tions  for  modeling  of  this  type  of  seismic  source.  Specifically,  we  can  simplify 
the  computations  implied  by  (B)  and  (9)  very  substantially  if  we  can  initially 
specify  the  rupture  velocity  function  in  a  manner  consistent  with  the  boundary 
conditions,  rather  than  be  forced  to  calculate  this  function  using  the  boundary 
conditions  (2)  and  (3)  jointly  with  the  field  relations  (B)  and  (9). 

Thus,  since  the  rupture  rate  can  be  simply  related,  by  equations  (14)  or 
(15),  to  the  intrinsic  medium  properties  and  the  dynamical  stresses  under  some 
rather  reasonable  general  assumptions,  the  following  approach  has  been  used 
in  modeling  earthquakes  for  purposes  of  predicting  the  properties  of  the  radi¬ 
ated  seismic  field  (see  also  Archambeau  and  Scales,  19B4,  for  additional 
details): 

(1)  An  initial  prestress  field  is  specified  and,  for  convenience  in  handling 
the  required  integral  evaluations  in  the  dynamical  solution,  the  field  is 
described  analytically  in  terms  of  the  associated  initial  displacement, 
expressed  in  vector  spherical  harmonics. 

(2)  A  failure  zone  geometry  and  time  history  of  evolution  is  specified  in 
terms  of  a  rupture  rate  function  Ur,  with  Ur  prescribed  in  a  form 
compatible  with  one  of  the  equations  (14).  (15)  or  (16). 

(3)  Static  solutions  for  the  equilibrium  displacement  field  of  an  inclusion, 
having  the  geometry  specified  in  step  (2),  are  obtained  either  numeri¬ 
cally  or,  when  possible,  analytically.  The  solution  is  expressed  in 
terms  of  a  vector  spherical  harmonic  expansion,  so  that  it  has  an 
analytical  expression  compatible  with  that  used  for  the  initial  dis¬ 
placement  field,  and  such  that  the  integrals  in  (B)  and  (9)  may  be  most 
easily  evaluated. 

(4)  The  expansion,  in  vector  harmonics,  of  the  rate  of  change  of  the 

equilibrium  field  (i.c.,  the  expansion  of  the  quantity  0u*(xo.fo)/  Ot  o 
where  Xo  and  t0  are  the  source  coordinates  and  time  variables)  is  then 
used,  in  the  integral  result  (6),  to  evaluate  the  initial  value  radiation 
field  This  field  is  expressed  in  vector  spherical  wave  functions, 

and  may  then  be  used  as  a  "starting  solution"  in  (9)  to  obtain  higher 
order  approximations  for  the  field,  given  by  In  all  of  this,  Ur  is 

employed  in  the  form  specified  in  step  (2). 

This  (approximate)  solution  procedure  can  be  specified  quantitatively 
using  the  formalism  developed  in  the  section  on  "Mathematical  Foundations". 
In  this  regard  the  static  initial  displacement,  as  well  as  the  equilibrium  dis¬ 
placement  field  in  the  vicinity  of  an  embedded  inclusion,  are  solutions  of  the 
static  elastic  equation 

V*u+  ^V(7u)s0 

where  o  is  Poisson's  ratio  for  the  material.  We  can  treat  the  medium,  as  an 
approximation,  as  uniform  and  extending  to  great  distances  from  an  inclusion 
of  variable  size,  such  that  its  dimensions  are  relatively  small  (or  very  small) 
when  computing  the  initial  displacement  field,  but  larger  and  of  variable  dimen¬ 
sions  when  computing  the  equilibrium  field  for  the  expanding  failure  zone.  We 
then  take  the  origin  of  coordinates  at  the  center  of  the  initial  (small)  inclusion, 
and  at  the  same  point  for  computation  of  the  growing  inclusion  solutions. 
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Further,  we  require  the  stress  field  to  assume  a  prescribed  functional  form  or 
value  at  great  distances  from  this  origin  such  that,  whatever  the  size  of  the 
inclusion,  the  equilibrium  stress  field  must  approach  this  value  (or  functional 
form)  at  great  distances  from  the  inclusion. 

In  this  case  all  the  equilibrium  displacement  fields,  including  the  initial  dis¬ 
placement  field  associated  with  the  prestress  state,  can  be  expressed  in  the 
form: 

*(*o)  =  [ato»H&»(*b)  +  fitmNH ,(r0)  +  7imFim(rb)]  +  u(0)(x) 

where  u^  is  the  displacement  associated  with  the  stress  prescribed  at  great 
distances  from  the  origin.  Here  Urn.  Him  &nd  are  exterior  vector  harmonics, 
which  are  eigenfunctions  of  the  Navier  equation  for  the  static  displacement 
field.  In  particular,  with  7  =  4(l  —  a). 

Hi,  vm+17cfc,(l>0.  *,) 


=  ^0l‘+8)  (l  (^0-  Fo)  “  (*  +  0  Pjm (^o-  ^o)j 


=  ro*  VriZTIT BfmK. Fo)  +  ifimfoo.  ?o) 

The  vector  functions  Pjm,  Htm  and  Can,  are  the  classical  vector  spherical  har¬ 
monics.  The  coefficients  a*n.  Ptm  and  Jim  are  constants  which  depend  on  the 
dimensions  of  the  inclusion  and  for  the  growing  inclusion,  therefore,  these 
coefficients  are  parametrically  dependent  on  time,  since  the  dimensions  of  the 
Inclusion  are  prescribed  by  the  rupture  rate,  Ur,  and  change  with  time  accord¬ 
ing  to  the  time  integral  of  Ur. 

The  initial  value  field  u*  is  defined  as  the  difference  between  the  equili¬ 
brium  field  for  the  growing  inclusion  (at  any  time)  and  the  initial  value  field. 
Thus,  using  the  previous  general  expression  for  these  equilibrium  fields  we 
have,  after  subtraction, 

«•<*)  =  2  [*AMfm(r)  +  AmNtai(r)  +  7*,F^(r)l  (18a) 

It*  ‘  1 


where 

a£,(*o)  =  Oim(fo)  - 

fiSnito)  =  fcm(*o)  -  Afi? 

7lm(*o)  =  7<m(*o)  ”7&\ 


(18b) 


These  coefficients,  as  noted,  are  (implicit)  functions  of  the  source  time  variable 
<0-  Here  etc.  are  the  constant  coefficients  associated  with  the  solution  for 
the  static  prestress  condition  prior  to  failure  zone  growth,  while  a*«(f 0)-  etc., 
are  the  coefficients  in  the  solution  for  the  inclusion  of  dimension  determined 
by  the  failure  zone  growth  rate  function  Ur. 

With  the  rupture  rate  prescribed,  or  approximated,  for  example,  by  a  func¬ 
tional  relation  such  as  (from  equation  16). 


where  v,  to  the  elastic  sheer  Telocity  and  |Ae|  is  the  magnitude  of  the  (partial) 
stress  drop  at  the  failure  boundary  while  |o,|  to  the  magnitude  of  ambient 
stress  at  the  same  point  just  before  failure,  then  we  can  prescribe  the  evolution 
of  failure  growth  and  consequently  the  time  dependence  of  the  coefficients  in 
(18).  (This  prescription  of  the  time  dependence  of  the  equilibrium  field  coeffi¬ 
cients  to  in  numerical  terms  when  the  failure  zone  inclusion  has  a  general 
geometrical  shape.)  Clearly,  if  we  take  the  stress  drop  magnitude  |&a|  to  be 
proportional  to  the  ambient  stress  |  a*  |  then  Uj »  is  a  constant  and  proportional 
to  the  shear  velocity.  This  approach  has  been  often  used  in  the  past,  and  some 
results  using  this  hypothesis  are  discussed  later.  On  the  other  hand  the  stress 
drop  can  (reaionebly)  be  taken  to  depend  non-linearly  on  the  ambient  stress, 
for  example  |Acr|  at  C0  |cr,  |  +  Cj  |  oB  |  *,  with  Co  and  Ct  constants.  In  this  case 
the  rupture  rate  would  vary  directly  with  the  ambient  stress  level,  that  is  we 
would  have 


V,  [c0  +  Cx  |o.|] 


hi  this  case  the  rupture  growth  would  speed  up  in  zones  of  high  initial  stress. 
As  will  be  pointed  out  later,  there  is  some  observational  evidence  that  this 
occurs. 

In  any  case,  based  on  plausible  hypotheses,  it  is  possible  to  define  a  rup¬ 
ture  rate  compatible  with  energy  and  momentum  conservation  and  to  thereby 
determine  a  time  history  of  the  failure  zone  development,  and  to  also  express 
this  time  history  in  terms  of  variations  for  the  equilibrium  field  expansion  coef¬ 
ficients  affa.  fli m  and  y 8*  in  the  equilibrium  field  u*. 

Now  we  have  that  the  dynamic  field  is  approximated,  to  first  order,  by  the 
relaxation  integral  term  given  in  (8).  hi  terms  of  the  Fourier  transformed  ver¬ 
sion  of  this  field,  using  the  results  and  inner  product  notation  given  in  the  sec¬ 
tion  on  "Mathematical  Foundations"  (equation  12),  we  have 


dirtf/1’  =  iuFt9  J  (pdtg vf.  C/)rJ 


where  /<§  denotes  Fourier  transform  operator  with  respect  to  t0,  the  source 
time.  Here,  using  the  previous  results  for  u*.  we  have  that 
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Farther,  the  transformed  Green’s  function  GK r,  r0:  «)  has  the  eigenfunction 
expansion  given  in  the  section  on  "Mathematical  Foundations"  (equation  9): 


„  f*7(r.nt>f)  f£(r0.  uf) 


f/(r.  nuf)  fT(r.  nuD 


with  f  and  pr  spheroidal  and  toroidal  vector  eigenfunction  for  a  layered 
spherical  earth  model,  expressible  in  terms  of  linear  combinations  of  the  classi¬ 
cal  vector  spherical  harmonics  P^n,  and  Ctm.  multiplied  by  Haskell- 
Thompson  matrices  involving  spherical  Bessel  functions. 

Using  these  expansions  in  the  approximate  solution,  in  (IB),  for  the 
dynamic  field  then  leads  to4 
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(20) 


Here  the  inner  product  form  is  used  to  denote  both  a  vector  inner  product  and 
a  functional  inner  product  corresponding  to  an  integration  over  the  volume 
K(f0);  that  is: 
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We  have  also  used  matrix  forms  in  (20)  to  be  able  to  write  the  results  in  more 
compact,  yet  explicit,  form.  Here  and  elsewhere  the  dagger  (t)  denotes  the 
transpose  of  a  matrix.  Thus,  in  (20),  the  expanded  matrix  and  inner  products 
when  written  out  give: 
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and  analogously  for  the  term  involving  the  toroidal  eigenfunctions. 

The  inner  product  integrals  over  V(fc)  can  be  evaluated  using  the  orthogo¬ 
nal  properties  of  the  vector  spherical  harmonics  P^,  Btm  and  CW  since  both 

Tht  summation  convention  it  applied  only  to  coordinate  indice*,  while  explicit  summation  is 
Specified,  when  appropriate,  for  the  "mode  indicee”,  l.m.n  and  l  \m\ 
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the  eigenfunctions  end  the  vector  solid  harmonics  are  linear  combinations  of 
these  functions.  In  particular,  for  a  thin  failure  zone  the  inner  products  can  be 
approximated  by  the  relation: 

<PHiW  6U6^-  (21) 

where  the  "weight  functions"  involve  only  integrals  over  the  radial  coordi¬ 
nate.  which  are,  in  general,  functions  of  the  source  time  variable  1 0  due  to  the 
time  dependence  of  the  region  of  integration  K(*o)-  Here  the  symbols  6u-  and 
a„u11-  are  Kroneclcer  delta  functions  having  zero  value  unless  1=1'  and  m  =  m.'. 
Thus  the  eigenfunction  expansion  of  the  radiation  field  for  a  growing  thin 
failure  zone  becomes,  using  (20)  and  the  previous  relation: 
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(22) 


The  structure  of  this  eigenfunction  expansion  for  the  source  radiation 
field  has  some  noteworthy  features,  in  particular  the  weight  functions 
and  r  are  all  indwpmdmt  of  the  source  properties  since  they  only  involve 
inner  product  integrals  of  the  spherical  earth  eigenfunctions  with  solid  vector 
harmonics.  (Thus,  for  a  given  earth  model,  they  only  need  be  calculated  once, 
and  cam  be  used  for  any  source  calculation.)  All  the  source  information  is  car¬ 
ried  in  the  coefficients  involving  time  derivatives  of  the  equilibrium  field  expan¬ 
sion  coefficients  a£,,  yffn.  These  coefficients  are,  therefore,  the  fundamen¬ 
tal  "descriptors'*  of  the  source  and,  in  an  inversion  procedure  designed  to  infer 
source  properties,  are  the  unknowns  to  be  determined. 

Ve  can  consider  the  entire  Fourier  transform  of  the  matrix  product  of  the 
weight  functions  times  these  coefficients  to  be  the  transform  of  a  scalar  func¬ 
tion.  That  is,  we  can  define 

[yr,&(fo)f 

a U(u)=f  *tf£n  sBtfftfto)  (23a) 
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and  rewrite  the  result  (21)  in  the  simple  form: 

aL^(o)  aL(») 


In  the  time  domain  the  inverse  Fourier  transform,  with  respect  to  the  time  vari¬ 
able  f,  is  easily  obtained  using  the  residue  theorem,  by  noting  that  the 
integrand  on  the  right  side  in  (22)  has  simple  poles  at  ±*of  and  ±nUi.  This 

gives: 

4nu(1,(r.t)  =  -  £  [tt&»(n«f)  f(rmwf)  cos* off 

+  a£m  (*«/)  fir^uD  COS* off]  (25) 

Clearly,  from  (24)  and  (25)  the  coefficients  (££,*  and  ££•*  are  the  quantities 
that  may  be  determined  from  observations.  Then  the  relations  in  (22)  are  those 
that  can  be  used  to  infer  first  order  source  properties  from  observations.  Here 
the  weight  functions  are  known  while  the  "equilibrium  coefficients"  are  to  be 
determined  in  such  an  ‘Inversion"  procedure.  On  the  other  hand,  computation 
of  the  "equilibrium  coefficients"  for  a  particular  source  geometry  and  time  evo¬ 
lution  is  all  that  is  necessary  to  completely  specify  a  theoretically  predicted 
dynamic  field  from  source. 

In  case  the  integrals  for  the  weight  functions,  defined  in  (21).  are  over  a 
region  which  is  fixed  in  time,  then  the  weight  functions  will  also  be  independent 
of  the  source  time  f0.  For  some  source  geometries  this  will  be  a  good  approxi¬ 
mation,  or  will  be,  in  other  special  cases,  exact.  In  any  case,  when  these  weight 
functions  are  independent  of  1 0.  or  nearly  so,  then  (23)  simplifies  to: 
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This  is  a  particularly  simple  result  which  is  easily  applied  for  the  inference  of 
source  properties  from  observations  or  for  the  direct  computation  of  theoreti¬ 
cally  predicted  radiation  fields  from  source  models. 

From  our  previous  considerations  of  the  elastic  wave  radiation  from  earth¬ 
quakes,  we  have  that  corrections  to  the  first  order  result  given  by  (24)  are 
obtained  by  computing  the  "scattering  integral",  given  in  equation  (12)  in  the 
section  on  "Mathematical  Foundations",  or  from  integrals  such  as  (9)  in  the 
present  section.  In  terms  of  the  inner  product  notation,  the  second  order  solu¬ 
tion  which  accounts  for  "scattering"  (and  absorption  of  energy)  on  the  failure 
surface  is  given  by: 

aJP(r.o)  *  ♦  Ft,  gH0ni)tV.  - 


Here  ell  the  functions  on  the  right  side  are  known  since  the  eigenfunction 
expansions  are  used  in  all  the  transformed  Creen's  tractions  and  displace¬ 
ments,  while  (23)-(25)  explicitly  specify  uj1*  and  from  their  definitions. 
While  this  form  can  be  manipulated  into  an  expansion  like  that  for  the  first 
order  field  from  the  source  (Archambeau  and  Stevens,  unpublished,  19B0) 
it  Is  not  worthwhile  to  give  these  complicated  results  here,  since  they  do  not 
need  to  be  used  in  the  subsequent  discussion.  That  is,  we  need  only  consider 
the  first  order  predictions  of  the  elastic  radiation  from  earthquakes  relative  to 
explosion  sources  in  order  to  quantify  their  differences. 

The  results  given  here  have  all  been  expressed  in  terms  of  expansions  in 
the  eigenfunctions  of  a  spherical  earth  model.  Similar  results  can  be  obtained 
for  plane  layered  earth  models,  wherein  the  expansions  are  in  terms  of  eigen¬ 
functions  in  plane  layered  geometries  (Archambeau  and  Stevens,  unpublished. 
10BO).  This  expansion  form  is  suitable  for  predictions  (and  source  property 
inversion  studies)  when  the  source  to  receiver  distances  are  less  than  a  few 
thousand  kilometers. 

.  Results  from  model  calculations  using  this  type  of  approximation  pro¬ 
cedure  have  provided  predictions  of  seismic  radiation  from  theoretical  sources 
that  are  reasonable  dynamical  approximations  of  earthquakes.  An  important 
feature  of  these  models  is  that  they  are  expressed  in  terms  of  the  fundamental 
physical  variable  controlling  the  dynamics  of  these  natural  sources  and  so  it  is 
possible  to  infer  both  the  basic  properties  of  the  seismic  field  and  its  variations 
in  terms  of  such  physical  parameters  as  those  describing  the  final  dimensions 
of  the  failure  zone,  its  rate  of  formation  and  prestress.  As  distinct  from 
kinematic  representations,  these  dynamical  solutions  do  not  involve  any 
assumptions  of  the  time  history  of  the  displacement  or  stress  field  at  the 
failure  surface,  or  'Yault  plane".  The  dynamical  modeling  does,  however, 
require  specification  of  the  initial  stress  field  and  some  general  assumptions 
regarding  the  physics  of  the  failure  process  involved,  such  as  the  assumption  of 
the  loss  of  shear  strength  in  the  material  after  failure. 

One  of  the  simplest  models  that  can  be  investigated  is  one  in  which  the  ini¬ 
tial  prestress  is  taken  to  be  uniform  in  the  medium  and  the  rupture  rate  to  be 
constant  and  near  the  medium  shear  velocity.  Various  rupture  zone  geometries 
have  been  considered  for  this  uniform  prestress  case,  all  of  them  corresponding 
to  some  type  of  ellipsoid  of  revolution  generated  by  simultaneous  expansion 
and  translation  of  spherical  failure  zones  (s .g.  Archambeau,  196B;  Minster, 
1973).  While  this  type  of  model  is  the  most  elementary  that  can  be  considered, 
it  does  allow  the  radiation  fields  to  be  evaluated  analytically,  so  that  the  depen¬ 
dence  of  the  seismic  field  on  the  various  source  parameters  can  be  evaluated 
and  described  by  analytic  means.  Most  important,  however,  is  the  fact  that  we 
can  expect  predictions  from  these  elementary  models  to  provide  a  description 
of  the  essential  first  order  features  of  the  seismic  radiation  associated  with 
natural  earthquakes.  Further,  we  also  expect  that  the  relaxation  source 
theory,  when  considered  in  the  limit  in  which  a  spherical  rupture  zone  is 
created  "supersonically",  will  proride  a  good  first  order  determination  of  the 
seismic  radiation  released  tectonically  by  an  underground  explosion  in  an  ini¬ 
tially  stressed  medium.  (For  details  of  the  modeling  of  explosion  induced  tec¬ 
tonic  release  see,  for  example,  Archambeau  and  Sammis,  1970;  Archambeau, 
1972.) 
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Analytical  Results  for  Simple  Relaxation  Theory  Models  of  Tectonic  Sources 

To  illustrate  the  procedure  used  to  construct  relaxation  theory  models  and 
to  provide  an  example  of  the  (analytical)  form  of  the  seismic  field  produced  by 
such  models,  we  will  summarize  the  mathematical  results  for  the  simple 
-  translating  and  expanding  spherical  failure  zone  models. 

i  In  order  to  simplify  matters,  while  nevertheless  maintaining  sufficient  gen- 
■[  erallty  for  purposes  of  describing  Important  properties  of  the  seismic  field,  we 
will  consider  a  homogeneous,  pur s  shear  initial  stress  field  as  a  representation 
of  the  (average)  prestressed  state  of  the  medium.  Further,  this  prestress  field 
will  explicitly  be  assumed  to  be  uniform  to  great  distances  from  the  coordinate 
origin  (is.,  the  point  of  initial  failure  in  the  medium)  and  so,  in  effect,  to  be  uni¬ 
form  to  infinity.  This  will  allow  us  to  extend  the  source  region,  and  specifically 
the  integration  over  the  zone  in  which  stress  relaxation  takes  place,  to  infinity 
and  to  thereby  simplify  the  Green's  function  representation  of  the  tectonic 
source.  (In  some  of  the  past  work  on  relaxation  source  representations,  for 
example,  by  Archambeau  (198B,  1972)  and  Minster  (1973),  the  region  in  which 
relaxation  of  a  homogeneous  prestress  was  allowed  to  contribute  to  the  seismic 
radiation  field  was  restricted  to  a  spherical  volume  of  radius  /&.  This  was  done 
on  the  grounds  that,  in  the  earth,  the  prestress  was  not  homogeneous  and  that 
the  zone  of  initial  stress  was  in  any  case  finite  and  not  infinite  as  is  implicitly 
assumed  when  a  homogeneous  initial  stress  is  used.  While  these  observations 
concerning  the  prestress  must  certainly  be  true,  the  approximation  of  the 
situation  in  the  earth  by  the  use  of  a  "cutoff'  in  the  homogeneous  prestress 
model  is  not  a  very  accurate  one.  In  particular,  subsequent  work  by  Stevens 
(1981)  has  shown  that  failure  induced  relaxation,  of  heterogeneous  prestress 
does  produce  seismic  radiation  having  spectra  that  are  similar  to  that 
predicted  using  a  homogeneous  prestress  along  with  the  "/i,  cutoff  factor",  but 
it  cannot  be  said  that  the  approximation  gives  an  accurate  representation  of 
the  low  frequency  radiation,  nor  does  it  show  the  high  frequency  complexities 
that  can  occur  when  the  prestress  is  strongly  heterogeneous.  Therefore,  in  the 
present  discussion  we  will  develop  the  radiation  field  predictions  for  the  com¬ 
pletely  homogeneous  prestress  case,  which  provides  a  good  first  order 
representation  for  earthquake  radiation  fields,  and  then  discuss  some  of 
Steven's  results,  among  others,  for  the  non-homogeneous  case,  contrasting 
them  with  results  from  the  homogeneous  case  in  order  to  show  how  a  spatially 
confined,  heterogeneous  initial  prestress  modifies  the  predictions.) 

For  the  case  of  a  homogeneous  pure  shear  prestress  the  mathematical 
description  can  be  developed  in  terms  of  the  dilatation  and  rotation  potentials, 
and  is  relatively  simple.*  Specifically,  the  dilatation,  and  vector  rotation 
fields,  x.  are  defined  by 

*The  un  of  theie  "physical  potentials"  is  limited  to  the  pure  shear  prestreas  ease  and  can¬ 
not  usually  he  used  to  represent  more  general  cases,  where  the  prestreas  is  arbitrary.  This 
is  because  these  potentials,  vhen  obtained  from  (87),  do  not  always  represent  the  entire 
elastic  field  is  the  static  case,  so  that  the  static  values  of  the  dilatation  and  rotation  poten¬ 
tials  do  not  fully  describe  prestress  or  equilibrium  stress  states  in  all  situations.  It  hap¬ 
pens,  in  the  pure  shear  prestress  case,  that  the  static  dilatation  and  rotation  potentials  can 
completely  represent  the  equilibrium  stress  state,  but  this  is  a  special  oaee.  In  the  general 
prestreas  case  the  displacement  field  itself  must  be  used  along  with  the  appropriate  dynam¬ 
ical  Croon's  function  integral  equation,  with  results  as  given  earlier.  If  potentials  are  to  be 
used,  however,  then  the  Lamd  potentials,  defined  by  U  *  7d  tTxf  are  required.  In  this 
latter  case  the  lamd  potentials  d  and  ^  satisfy  ordinary  wave  equations  and  the  analysis 
proceeds  in  the  same  way  as  it  does  for  the  dilatation  and  rotation,  with  formal  results  that 
are  only  different  by  constant  factors. 


Xs\Xi>X*'Xa)s  jVxu 
*»  =  V-u 


The  equation  of  motions  for  an  isotropic  and  homogeneous  medium,  with  negli* 
gible  body  forces  acting,  is 

9fu  *  u/V*4  -  2u/(V  x  x)  (2Ba) 

with  Vp  and  v,  the  compressional  and  shear  velocities  in  the  medium.  There¬ 
fore  the  equation  of  motion  can  be  used  to  compute  the  acceleration  from 
these  potentials.  If  the  Fourier  transform  with  respect  to  time  is  applied  to  this 
relation,  then  the  relation  is: 


(28b) 


where  kp  *  o/Vp  and  k,  •  u/va  are  wave  numbers  for  compressional  and  shear 
waves.  Thus  displacement  transforms  can  be  computed  from  the  potential 
transforms  using  this  relation. 

It  is  not  difficult  to  show  that  the  equation  of  motion  is  satisfied  if  the 
Cart  Brian  components  of  x  and  the  scalar  potential  all  satisfy  wave  equations 
of  the  same  form  (Archambeau,  1988).  In  particular,  the  potentials  must 
satisfy: 


where  a  *=  1,  2,  3,  4  and 


V*X*  ~  7T  =  0 


(Xa)  =  (Xi.  Xa-  Xs<  X*) 


is  used  to  denote  any  one  of  the  four  potentials.  Here  also 

(v.)  s  (v„i/„  v„vp) 

is  the  corresponding  shear  or  compressional  velocity  associated  with  the  indivi¬ 
dual  potentials. 

The  advantage  of  the  use  of  these  potentials  is  that  the  eigenfunction  solu¬ 
tions  of  the  wave  equations  in  (29)  ere  scalar  wave  functions  with  well  known 
properties  (see,  for  example,  Morse  and  Feshbach,  1953)  and  are  easily  manipu¬ 
lated,  as  opposed  to  the  much  more  complicated  and  cumbersome  vector  wave 
functions  associated  with  the  elastic  equations  of  motion.  Further,  as  will 
become  apparent,  the  radiation  field  is  conveniently  separated  into  purely 
compression  waves  (x4)  and  shear  waves  (x)  by  this  potential  decomposition. 

The  Green's  function  solution  of  (29)  corresponding  to  the  volume  relaxa¬ 
tion  source  term,  or  initial  value  integral  term  denoted  by  in  equation  (8), 
is  simply  (Archambeau,  1988): 

1  *r  ..  r  8I<*>  . 
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where  I*BI(r,f;ib.<o)  !■  the  scalar  Green’s  function  for  the  wave  equation.  The 
volume  integral  is  over  V(t0),  the  volume  outside  the  failure  zone.  This  result  is 
completely  analogous  to  the  first  order  solution  (u^1))  discussed  earlier  and 
given  in  equation  (8).  The  potential  xf  appearing  in  this  integral  solution 
corresponds  to  the  change  In  the  equilibrium  value  of  the  potential  x,-  as  a 
function  of  the  source  time  variable  to<  due  to  the  (spontaneous)  creation  and 
growth  of  the  failure  zone.  This  function  is  completely  analogous  to  the  equili¬ 
brium  displacement  function  a*,  and  can,  in  the  present  case,  be  computed 
from  it,  using  the  potential  definitions/ 

The  limit  on  the  integral  over  the  source  time,  is  t*  =  t  +  e,  with  t  >  0 
and  small,  where  e  is  introduced  simply  to  avoid  singular  points  in  the  Green’s 
function.  However,  if  1 0  >  r0,  with  r0  the  total  time  duration  for  formation  of 
the  failure  zone,  then  x£  does  not  vary  for  such  values,  since  the  failure  zone 
does  not  change  or  grow,  and  then  the  derivative  of  x£  will  vanish.  Therefore 
the  integral  over  the  source  time  1 0  has  an  upper  limit  equal  to  r0  when  t  >  r0. 
In  general  we  will  consider  the  case  in  which  the  observer  time  t ,  measured 
relative  to  the  beginning  of  failure,  is  larger  than  the  time  To*  In  this  case  we 
will  replace  the  limit  by  Tq  and  explicitly  note  that  this  representation  is  valid 
for  t  >  Tq. 

The  use  of  the  index  (a)  with  the  Green’s  function  is  to  indicate  that 
slightly  different  Green's  functions  must  be  used  for  the  rotation  potential 
components,  i  sl>  8,  3.  than  is  used  for  the  dilatation  potential  In  par¬ 
ticular,  to  compute  the  direct  radiation  field  from  a  relaxation  source,  we  can 
use  the  infinite  space  Green's  function,  and  in  this  case: 


*0'  *o) 


r* 


where  6(x )  denotes  a  Dirac  delta  function  and 

r*  =  |  r  -  r„  | 


t*  =  t  -t0 

with  r  and  t  receiver  coordinates  and  time.  Here  the  Green's  function  for 
a  *  1,  2,  3  involves  the  shear  velocity  in  r*B\  while  the  compressional  velocity  is 

used  when  a  *  4. 

The  Fourier  transformed  integral  solution  corresponding  to  (30)  is 
obtained  by  taking  transforms  with  respect  to  the  observers  time  t .  Thus 

=  f  e*«  dt 


Is  the  transformed  potential.  Applying  this  operation  to  (30),  and  taking 
account  of  the  properties  of  the  Dirac  delta  function,  we  have  for  t  >  r0: 
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dr0;  t  >  r0 


(31) 


is  noted  earlier,  the  dilatation  and  rotation  potentials  do  not  always  completely  describe 
the  equilibrium  state  a*  when  the  potential  definitions  in  (Z?)  are  used  for  state  field  com¬ 
putations. 


Ibis  spectral  representation  is  most  useful  for  the  evaluation  of  the  potentials, 
as  compared  to  the  time  domain  representations,  and  (31)  is  generally  the  form 
that  has  been  used.  (See,  for  example,  Archambeau,  1964,  1968  and  Minster, 
1973.) 

The  integral  representation  of  the  potentials  can  be  evaluated  using  a  pro* 
eedure  which  is  equivalent  to  that  described  earlier  for  approximate  solutions 
for  the  source  displacement  field.  This  procedure  is,  as  before,  an  approxima¬ 
tion  to  the  complete  dynamical  solution,  wherein  we  specify  the  rupture 
growth,  rather  than  solve  for  it  using  the  boundary  conditions  and  the  equation 
jgf  state  for  the  material.  (Further,  for  this  example,  only  the  initial  vqjue  field 
jgi1)  will  be  obtained,  with  the  scattered  field  contributions  given  by  x£,n)  with 
ns2  completely  neglected.) 

One  way  of  specifying  failure  sones  is  indicated  in  Figure  1.  In  particular, 
the  simple  spherical  zone  can  be  made,  to  expand  and  translate  in  such  a  way 
that  the  time  evolution  of  a  natural  or  explosion  induced  failure  process  can  be 
approximated.  (Clearly,  a  thin  disk  or  ellipsoid  would  be  a  more  preferable 
geometric  shape  to  use  for  simulation  of  an  earthquake  failure  zone,  rather 
than  a  sphere.  However,  the  moving  expanding  sphere  has  the  advantage  of  giv¬ 
ing  rjelstively  simple  closed  form  analytical  results  for  the  potentials,  while  the 
other  geometries  produce  much  more  complex  and  cumbersome  results. 
Further,  as  indicated  hy  Minster  (1973),  when  comparisons  are  made,  the  radia¬ 
tion  field  predictions  for  models  with  "thin  failure  zone"  geometries  do  not 
differ  in  any  fundamental  way  from  those  using  a  tangentially  expanding 
sphere  with  the  same  final  dimensions  and  growth  characteristics.) 

-  t  For  a  spherical  zone  at  any  time  f0,  we  can  parameterize  the  failure  zone 
in  terms  of  two  functions  of  time:  d(t0).  a  function  describing  the  movement  of 
the  center  of  the  sphere  relative  to  a  fixed  coordinate  system,  and  tf(f0),  the 
time  dependent  radius  of  the  sphere,  as  indicated  in  Figure  1.  By  adjusting  the 
time  dependence  of  these  two  functions,  we  can  obtain  time  evolutions  such  as 
those  illustrated  in  the  three  lower  insets  in  Figure  1.  (This  parameterization  is 
equivalent  to  a  full  specification  of  the  rupture  velocity  function,  Ur,  as  a  func¬ 
tion  of  time  and  position.)  Of  the  three  models  shown,  the  uniformly  expanding 
sphere  has  been  used  to  model  tectonic  release  effects  associated  with  explo¬ 
sions  in  prestressed  media  (e .g.,  Archambeau,  1972)  while  the  tangentially 
expanding  sphere  and  moving  sphere  models  have  been  used  to  approximate 
earthquake  failure  zones  by  Minster  (1973)  and  Archambeau  (19BB),  respec¬ 
tively.  Of  the  latter  two  model  types,  the  tangentially  expanding  model  leads  to 
the  most  accurate  approximate  representation  for  earthquakes. 

Now  the  equilibrium  displacement  field  In  the  elastic  medium  surrounding 
a  spherical  cavity  in  a  homogeneously  prestressed  medium  is,  from  Landau  and 
Lifshitz  (1959)  for  example: 
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Here,  is  the  constant  prestress  existing  prior  to  the  introduction  of  the  cav¬ 
ity,  o  is  Poissons  ratio,  p,  is  the  elastic  rigidity  of  the  medium  and  /?0  is  the 
radius  of  the  spherical  cavity.  If  we  specify  the  prestress  to  be  a  constant  pure 
shear  field,  with: 
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figure  1.  Coordinate  eyatems  uaad  to  deacrlbe  th*  geometry  of  a  tranalatlng 

end  expanding  aphere.  The  primed  coordinate  ayatem  moves  with  the  center 
of  Che  aphere  while  the  fixed  coordinates, with  origin  (0),are  used 
to  locate  aource  polnta  Q  In  the  system  r  ,  8  ,  *  and  observer  points 
f  in  the  system  r,  8,  8.  The  lower  three  insets  illustrate  three 
geometric  models  for  failure  tones  using  translating/expanding  spheres. 
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then  this  solution  can  be  used  to  represent  the  equilibrium  displacement  field 
around  a  spherical  inclusion  in  which  the  shear  strength  vanishes  (jx  *  0  inside) 
and  in  which  the  compressibility  matches  that  of  the  medium  outside  the  inclu¬ 
sion.  In  this  case,  the  associated  potentials  xS  can  be  shown  from  (27)  to  have 
the  form,  in  the  coordinates  r'o,  d’o.  p'o  in  the  frame  O'  in  Figure  1: 
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where  the  coefficients  and  a  a  1,  2.  3,  4.  are: 
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with  d««  a  o  if  a  ^  4.  and  dM«  »  1  if  a  a  4.  Here  the  potential  function  index  (a) 
numbers  the  rows  vertically  from  1  to  4  in  the  matrices,  while  m  =  0,  1,  2  is  the 
column  index  running  from  left  to  right.  (This  result  is  given  by  Archambeau 
(1984,  1968)  with  a  multiplicative  factor  of  1/2  missing  from  the  rotation  poten¬ 
tial  coefficients.  Minster  (1973)  gives  the  correct  results,  but  omits  a  minus 
sign  from  the  coefficient  a^.  The  result  given  here  is  from  Harkrider  (per¬ 
sonal  communication),  who  gives  the  general  result  for  erf  ft  p  o$  P  oj§)  P  0  In 
this  general  case  the  following  additional  coefficients  are  non-zero  in  the 

matrices  a$  and  6$:  o|?  a  j  (aj?>  -  o$>),  a  (2cJ§>  -  cr|?>  -  0$), 

•It1  *  J  (gR1  “  qJP).  6il)  =-“  \  (0$  -  0$$).  ®nd  a  (<$>  -  q|?>).  In  the 
case  of  a  pure  shear  initial  stress,  these  coefficients  all  vanish.) 

With  this  description  of  a  "spherical  inclusion  equilibrium  field,  we  can 
express  the  equilibrium  field  at  any  source  time  f0>  for  *ny  of  the  failure  zone 
evolutions  in  Figure  1,  by  using  an  appropriate  radius  function  R(t 0).  The  field 
so  expressed  will,  however,  be  in  the  moving  coordinate  system,  with  origin  O' 
and  spherical  coordinates  r\  d',  p*. 


The  movement  of  this  coordinate  frame  is  a  simple  translation  relative  to 
the  fixed  coordinates  with  origin  0,  and  is  described  by  the  function  d(f0)- 
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Thus,  to  use  this  description  of  the  equilibrium  field  x»  the  integral  represen¬ 
tation  of  the  dynamic  radiation  field,  in  (31).  it  is  necessary  to  either  express 
the  field  Xm  terms  of  the  coordinates  of  the  fixed  reference  frame  (r0.  do.  po) 
and  evaluate  (31)  as  it  stands,  or  to  evaluate  the  integral  representation  for  the 
dynamic  field  in  the  moving  frame,  with  source  coordinates  (f*Q ,  do*  Po)  and 
observer  coordinates  (r',  d',  p'),  and  to  then  transform  the  dynamical  results 
back  to  the  fixed  reference  frame.  The  former  procedure  was  used  by  Archam- 
beau  (1968),  while  Minster  (1973)  used  the  latter.  We  shall  follow  Minster’s  pro¬ 
cedure  since  the  results  are  somewhat  more  general  and  compact  than  the 
alternate  expansion. 

The  first  step  in  evaluating  (31)  is  to  express  the  infinite  space  Green’s 
function  in  the  integrand  in  terms  of  spherical  wave-functions  referenced  to 
the  coordinate  system  which  moves  with  the  expanding  spherical  zone  of 
failure.  This  is  a  classical  result  (e.y.,  Morse  and  Feshbach,  1953),  and  has  the 
form: 
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where  (r‘,  d\  p')  are  source  point  coordinates  in  the  moving  frame,  correspond¬ 
ing  to  the  point  Q  in  the  figure.  The  angle  7'  is  measured  between  the  vectors  r* 
and  i £.  (The  vector  t,  from  the  origin  O'  to  the  observer’s  point  P  is  not  shown 
in  Figure  1,  but  vectors  r  and  i*  to  the  point  P  have  the  same  relationship  to  the 
coordinate  systems  at  0  and  O'  as  do  the  vectors  r0  and  ro  for  the  source  point 
Q).  The  Legendre  function,  involving  the  angle  y\  has  the  expansion  (Morse  and 
Feshbach,  1953) 

ft(cos7‘)  *  £  (2  -  dmo)  iT(c osd  )  PT( cosdi)  cosm(p’  -  pi) 


which  involves  the  basic  observer  and  source  point  coordinate  angles  d'.  p'  and 
pi-  Here  denotes  a  Kronecker  delta  function.  Hence,  these  two  expan¬ 
sions  allow  us  to  express  the  Green's  function  in  separated  form  in  terms  of  the 
source  and  observer  coordinates  in  the  moving  frame.  Now,  using  thij  Green’s 
function  expansion  and  the  solution  for  Xa>  the  equilibrium  potential  field  for 
the  spherical  inclusion  as  expressed  in  (32),  we  get,  from  (31),  the  integral  solu¬ 
tion: 
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tor  the  dynamic  field  due  to  relaxation  effects,  for  observer  times  t  >  To.  where 
To  is  the  time  necessary  for  the  complete  formation  of  the  failure  zone.  Here 


we  have  factored  out  the  time  dependence  term  from  the  multiple  coefficients 
and  b$  in  (32).  ao  that  ajjg}  and  b^Xj}  are  the  time  independent  coefficients 
defined  by: 

(afcidc)]  •  *(«.)  («ft>) 

Observing  that  the  orthogonality  of  the  Legendre  functions  is  such  that: 
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The  spatial  integrals  involving  the  spherical  Bessel  function  and  Hankel  func¬ 
tion  are  standard  integrals  which  may  be  evaluated  from: 

/ A(*. r„)  rjir.  =  (t<oy_V  -  Jj^TT 
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Using  these  relations  and  the  Wronskian  relation  for  the  spherical  Hankel  and 
Bessel  functions  given  by  (e.p.,  Abramowitz,  1964): 
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then  we  can  put  the  result  for  Xe1^  into  the  form: 
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where,  in  the  final  integral  term,  we  have  reintroduced  the  original  multiple 
coefficients,  afiH(to)  end  bjgifct 8),  associated  with  the  equilibrium  field  Xa- 

Now  we  observe  that  the  final  integral  is  not  a  propagating  field  but  actu¬ 
ally  corresponds  to  the  Fourier  transform  of  the  changing  initial  value  field  xS • 
times  the  Fourier  transform  of  a  step  function  which  accounts  for  the  factor  of 

multiplying  the  integral.  Thus  this  term  results  from  the  time  dependence 

of  the  initial  value  field  as  the  failure  tone  is  created  and  corresponds  to  the 
change  in  the  potential  due  to  the  changing  reference  state.  These  reference 
state  changes  are,  by  choice,  measured  relative  to  the  (fixed  or  static)  final 
equilibrium  state  in  the  formulation  of  the  original  integral  solution  given  in 
(22),  and  so  this  term  accounts  for  the  fact  that  this  final  equilibrium  is  not 
achieved  until  dynamic  relaxation  of  the  stress  has  taken  place  everywhere  in 
the  medium.  Hence,  the  term  functions  to  keep  track  of  the  instantaneous 
equilibrium  state  relative  to  this  final  reference  state. 

The  first  integral  term  in  the  expansion  for  xP  *  propagating  wave  field, 
in  view  of  the  presence  of  the  spherical  Hankel  function,  and  corresponds  to 
the  radiation  field  that  would  be  measured.  Ve  are  therefore  interested  in  the 
detailed  form  of  this  field  from  the  dynamical  standpoint,  rather  than  in  the 
second  integral  term,  and  will  simply  omit  the  second  integral  from  further  con¬ 
sideration.  ^ 

Use  form  of  the  dynamical  term  in  is  that  of  a  quadrupole  radiator  in  a 
moving  coordinate  frame  (r',  p').  However,  the  coordinates  in  the  moving 

frame  implicitly  depend  on  the  source  time  t0,  since  the  moving  frame  is 
defined  to  have  its  origin  at  the  center  of  the  expanding/translating  sphere 
with  radius  and  position  depending  on  f0  111  *ome  parametric  fashion.  To  bring 
out  this  dependence  explicitly,  and  to  evaluate  the  integral  over  the  source 
time  t0,  we  need  to  express  the  wave  field  in  the  fixed  coordinate  system 
(r,  p)  indicated  in  Figure  1.  That  is,  we  must  transform  the  expression  for 
the  quadrupole  radiation  field  in  the  moving  frame  to  an  equivalent  multipole 
field  in  a  fixed  reference  frame.  To  do  this  we  can  use  the  addition  theorem  for 
spherical  wave  functions,  as  employed  by  Minster  (1973).  In  particular,  for  the 
translation  by  d(t9)  of  the  moving  frame  along  the  s  axis  of  the  fixed  reference 
frame,  as  shown  in  Figure  1,  Minster  shows  that  the  spherical  wave  function  of 
order  1=2  transforms  between  the  moving  and  fixed  frames  according  to: 
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Here  the  coefficients  ct(i/,f|2,  m)  are  related  to  Cle bach- Gordan  coeffi¬ 
cients.  where  in  general 

Cj(i /,  f  |n.  m)  =  +  1)(2Z  ♦  1)* (2n  -4- 1)"*  •  (l  um  0  |  n  m)  (i  i/OO  |  n  0) 

The  symbols  (lvmO|nm),  etc.,  on  the  right  are  the  Ciebsch-Gordan  coeffi¬ 
cients,  and  since  they  are  tabulated  (e .g.,  see  Edmonds  (1957)  or  Abramowitz 
and  Stegun  (1965)),  their  computation  presents  no  problem.  The  relation,  as 
noted,  is  valid  when  r'  >  d(t0).  When  r'  <  d(t0).  corresponding  to  the  very  near 
field  distance  range,  then  another  form  of  this  expansion  must  be  used  (see 
Minster,  1973  for  details).  However,  our  main  interest  is  in  the  far,  or  inter¬ 
mediate,  distance  range,  so  that  the  expansion  given  here  is  sufficient. 

Substitution  of  the  expansion  for  the  “moving  quadrupole  field"  in  terms  of 
a  fixed  reference  frame,  wherein  a  translational  aspect  of  rupture  zone  growth 
is  included,  gives  finally: 
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We  have,  therefore,  that  the  radiation  field  due  to  the  relaxation  effects 
associated  with  failure,  when  expressed  in  terms  of  the  potentials  xi!) 
corresponding  to  the  scalar  dilatation  (a  =  4)  and  vector  rotation  (a  =  1,  2.  3), 
is  a  multipole  field  whose  coefficients  and  B&^  are  (independent)  functions 
of  frequency,  which  depend  on  the  growth  rate  of  the  failure  zone  through  the 
functions  J?(i0)  and  d(t0 )  for  the  particular  parameterization  used  here,  as  well 
as  upon  the  initial  stress  state  of  the  medium  and  its  elastic  properties. 
Clearly,  the  case  considered  here  is  special,  since  it  has  been  assumed  from  the 
onset  that  the  prestress  was  uniform  and  pure  shear  and  that  the  failure  zone 
was  of  a  type  described  by  a  translating  and  expanding  sphere.  However,  it  is 
nevertheless  true  that  the  radiation  field  from  any  source  can  be  expressed  in 
the  form  of  (34a),  and  it  is  only  the  multipole  coefficients  that  change  form  as  a 
function  of  the  detailed  nature  of  the  source. 

The  specific  form  of  the  result  obtained  in  (34)  will  depend  on  the  choice  of 
the  functions  R(t0)  and  d(f0).  Once  this  choice  has  been  made,  the  evolution  of 
the  failure  zone  is  completely  defined  and  the  integral  in  (34)  can  be  evaluated. 
This  integral  is  seen  to  be  a  Fourier  transform,  involving  products  of  spherical 


Bessel  functions.  In  general  it  is  evaluated  numerically  by  the  Fast  Fourier 
Transform  (FFT)  method.  (See  Minster,  1973,  however,  for  methods  of  analytical 
evaluation.)  Limiting  analytical  forms  of  the  integral,  at  large  or  small  values  of 
frequency  v  for  example,  can  be  obtained  without  much  difficulty  in  specific 
cases  of  interest  and  Minster  (1973)  provides  examples.  We  shall  use  these 
results  later,  in  the  section  on  "Scaling  Laws”,  to  provide  a  description  of  the 
seismic  radiation  in  terms  of  basic  event  parameters,  such  as  rupture  rate, 
prestress  and  failure  zone  dimensions. 

The  case  of  rupture  zone  evolution  illustrated  in  figure  1  are  among  those 
that  have  been  used  to  model  tectonic  radiation  effects  from  both  earthquakes 
fluid  underground  explosions.  As  mentioned  earlier,  of  the  models  illustrated  in 
figure  1,  the  tangentially  growing  spherical  failure  zone  model  produces  radia¬ 
tion  field  predictions  that  are  considered  to  most  accurately  approximate 
earthquake  radiation  fields.  In  this  case  the  choice  for  0)  fluid  d(tD)  is  con¬ 
strained  by: 

d(t0)  *=  R(t0) 

andi  for  uniform  growth,  consistent  with  the  uniform  prestress,  then 

d(*o)  +  *(«d)  *  2f?(f0)  =  Lfcf0 

is  the  necessary  choice  based  on  our  earlier  considerations  of  proportionality 
of  rupture  rate  to  prestress  level,  while  Ur  is  a  constant  rupture  rate.  In  this 
case  the  multipole  coefficients  for  the  potentials  in  (34)  are  given  by: 
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where  Jt0  i*  the  final  radius  of  the  spherical  failure  zone.  It  is  important  to 
notice  that  the  rupture  rate  appears  in  the  result  through  the  ratio  (va/  Ur) 
and  so  the  radiation  field  is  only  affected  through  the  ratios  of  the  intrinsic 
elastic  velocities  to  the  rupture  velocity.  This  also  shows  that  P  and  S  waves  will 
have  differently  shaped  radiation  spectra,  since  this  ratio  and  the  integral  limit 
have  different  values  for  the  two  wave  types. 

The  integral  in  (35)  can  be  evaluated  as  a  Fourier  transform,  or  since  the 
Bessel  functions  are  a  finite  series  of  sines  and  cosines  multiplied  by  polynomi- 
sds  in  (l/z),  the  integrand  can  be  expanded  to  a  finite  series  and  evaluated 
term  by  term.  The  results  will  be  discussed  in  later  sections. 

The  case  of  a  uniformly  expanding  spherical  failure  zone  illustrated  in  Fig¬ 
ure  1  is  also  of  considerable  interest,  not  only  because  the  result  turns  out  to 
be  simple  and  easy  to  understand,  but  also  because  it  can  be  used  to  model  tec¬ 
tonic  release  effects  associated  with  an  explosion.  In  this  case,  d(*o)  =  0.  end 
the  results  in  (34)  reduce  to  (see  also  Archambeau,  1972): 
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with  i?0  the  final  failure  zone  radius,  B(t0)  =  Ur  to  and  Ur  <  vm. 

If  the  rupture  rate  is  less  than  the  elastic  velocities  in  the  material,  then 
these  results  apply.  If  the  failure  process  is  driven  by  a  shock  wave  from  an 
explosion  however,  then  the  rupture  rate  can  be  greater  than  the  intrinsic  elas¬ 
tic  velocities  and  then  ws  must  use  the  "instantaneous  source"  integral 
representation,  equivalent  to  equation  (6)  discussed  earlier.  In  this  case  we  get 
the  somewhat  different  result: 


(38c) 


where  the  rupture  rate  Ur  vm  is  assumed  constant  and  Bo  is  the  radius  of  the 
final  failure  zone.  (This  result  follows  from  (36b)  if  we  insert 
J?(t0)  =  BoH(tc  -  Bo/  *ith  H  denoting  a  step  function.)  These  results  can 
also  be  expressed  in  terms  of  elementary  functions  if  we  use  the  identify: 
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In  any  cue,  the  results  in  (36)  show  that  the  seismic  radiation  to  be  expected 
from  the  production  of  a  spherical  failure  zone  in  a  uniformly  prestressed 
medium,  with  the  prestress  being  a  pure  shear  field,  is  a  pure  quadrupole  field, 
having  a  simple  double  couple  force  equivalent.  The  characteristics  of  the  radi¬ 
ated  spectrum  are  also  very  simple,  and  in  later  sections  we  shall  illustrate  the 
important  properties  of  these  spectra  using  these  and  similar  results. 

Results  and  Predictions  from  Relaxation  Theory  Modeling 

Figure  2  shows  results  obtained  by  Minster  (1073)  for  the  displacement 
spectra,  for  compression  (P)  and  shear  waves  (SV  and  SH),  from  a  simple,  uni¬ 
form  prestress,  relaxation  theory  model  employing  a  tangentially  expanding 
spherical  failure  zone.  These  results  correspond  to  the  initial  value  field,  and 
do  not  include  the  effects  of  scattering  from  the  failure  surface. 

The  important  first  order  characteristics  of  the  spectra  illustrated  here 
•re  however  (l).  the  large  difference  in  P and  Suave  spectral  levels,  with  the  S 
wave  production  nearly  a  order  of  magnitude  larger  at  all  frequencies.  (2),  the 
difference  in  'bomer  frequencies"  (that  is,  the  frequency  at  which  the  spectral 
roll-off  at  high  frequencies  clearly  begins )  between  the  P  and  S  wave  spectra. 
with  the  S  wave  spectra  having  a  significantly  lower  comer  frequency  (The 
spectral  holes,  evident  in  the  spectra  at  one  of  the  azimuths,  are  primarily  due 
to  neglect  of  the  boundary  scattering  terms  in  the  solution  and  are  essentially 
removed  when  scattering  is  taken  into  account  in  the  solution.)  (3),  the  far 
field  spectra,  for  all  wave  types,  tends  to  flatten  at  frequencies  lower  than  the 
comer  frequency  far  each  wave  type  (P  or  S),  this  occurring  when  the  prestress 
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ia  uniform  (or  marly  uniform)  to  large  distances  from  the  failure  none.  These 
three  properties  of  the  radiation  spectra  are  important  distinctive  characteris¬ 
tics  of  tectonic  relaxation  sources  and  are  among  those  most  clearly  observed 
for  earthquakes. 

figure  3  illustrates  the  effects  of  rupture  rate  on  the  shape  of  the 
congressional  (P)  and  shear  (SH)  wave  spectra  from  a  model  similar  to  that  of 
Figure  1.  Here,  in  addition  to  the  previous  three  observations  concerning  the 
wave  spectra,  it  is  clear  that  the  corner  frequencies  and  high  frequency  spec¬ 
tral  slopes,  for  both  P  and  S  waves,  are  strongly  affected  by  changes  in  the 
value  of  the  rupture  rate,  denoted  as  Vs  in  the  figure.  In  particular,  this  figure 
illustrates  two  additional  important  features,  namely:  (4)  The  comer  frequencies 
far  both  P  and  S  waves  are  directly  proportional  to  the  rupture  rate.  Further, 
(5):  The  high  frequency  spectrum  of  the  Pwave  varies  asymptotically  as  u-9  for 
high  rupture  velocities  near  the  shear  velocity  of  the  medium,  while  the  S  wave 
spectrum  falls  off  as  at  high  frequencies  for  such  high  rupture  rates,  fbr 
lower  rupture  velocities,  both  the  P  and  S  wave  spectra  show  quite  wide  fre¬ 
quency  bands,  immediately  above  their  corner  frequencies,  where  the  spectra 
vary  as  tt~x  to  u~*  before  assuming  steeper  slopes  at  yet  higher  frequencies. 

Comparisons  between  spectra  with  different  failure  zone  dimensions,  such 
as  those  of  Figures  (2)  and  (3),  shows  that:  (6)  The  spectral  levels  all  scale  as 
the  cube  of  the  failure  sons  length  when  the  failure  sane  surface  area  increases 
as  the  square  of  the  length  (which  is  usually  the  case  far  small  events).  The 
spectral  levels  increase  as  the  square  of  the  failure  sane  length  when  the  length 
is  allowed  to  increase  without  changes  in  the  other  failure  sane  dimensions 
(which  is  most  likely  for  very  large  earthquakes).  In  addition:  (7)  The  comer 
frequencies  of  the  P  and  S  wave  spectra  vary  inversely  with  the  (final)  length  of 
the  failure  cone. 

Figure  4  ahows  compressional  and  shear  wave  radiation  patterns  displayed 
at  three  different  frequencies.  The  important  radiation  field  properties  illus¬ 
trated  here,  which  can  be  added  to  those  previously  described,  are:  (B)  The 
radiation  patterns  at  low  frequencies,  that  is,  below  the  comer  frequencies  for 
each  wave  type,  are  dominated  by  a  quadrupole  term  (having  a  double  couple 
force  equivalent),  while  at  frequencies  above  the  corner  frequency,  for  each 
wave  type,  the  patterns  are  strongly  distorted  by  higher  order  multipole  contri¬ 
butions  which  are  associated  with  rupture  propagation.  The  distortion  of  the 
patterns  is  such  that  the  patterns  rotate  toward  the  direction  of  rupture  propa¬ 
gation  of  frequencies  above  the  "comer  frequency "  and  more  of  the  high  fre¬ 
quency  energy  is  radiated  inside  a  cone  along  the  rupture  plane  in  the  direction 
of  rupture  propagation. 

The  properties  described  constitute  the  most  basic  and  robust  features  of 
the  seismic  radiation  predicted  for  earthquakes,  and  they  have  a  reasonably 
good  correlation  with  observations.  In  cases  where  we  can’t  be  absolutely  con¬ 
fident  tbat  the  earthquake  data  is  consistent  with  a  predicted  feature,  it  is 
because  other  effects,  such  as  anelastic  attenuation,  scattering,  or  ordinary 
wave  propagation  effects  in  an  uncertain  structure,  preclude  definitive  verifica¬ 
tion  due  to  ambiguity  of  interpretation,  or  simply  because  the  observed  data  is 
too  limited  or  is  too  contaminated  by  noise. 

In  addition  to  the  characteristics  so  far  described,  there  are  effects  that 
will  arise  from  strong  departures  from  uniform  prestress  conditions.  This  will 
•how  up,  as  was  already  indicated,  in  variations  in  the  rupture  velocity  and 
•tress  drop  during  the  growth  of  the  failure  zone.  If  these  variations  are  strong 
enough,  then  we  can  simply  view  the  event  as  a  superposition  of  several  nearly 
constant  stress  drop  events  each  with  a  rupture  rate  scaled  appropriately,  as 
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Figure  3.  Theoretical  displacement  spectre  from  a  relaxation  source 
model  of  a  1  km  length  earthquake  shoving  effects  of 
variations  in  the  rupture  rate  on  the  spectral  shapes  for 
f  and  S  waves.  From  Minster  (1973). 
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Coaprasslonal  (?)  and  shear  v.v.  CSV* nd  SH)  radiation  pattern*  from • 
relaxation  thaory  sod cl  of  *n  earthquake  shoving  the  quadrupola  (N  •  2) 
patterns  and  the  sums  of  quadrupola  plua  higher  order  ■ultipole  tarns  (up  to 
V  -  4  and  K  -  6)  at  three  different  frequencies.  The  patterns  are 
vertical  sections,  with  T  denoting  tha  "take-off"  angle  from  the  *furce 
hypocenter  and  with  the  rCpture  velocity  vector  indicated  by  VR.  The 
source  nodal  has  a  "corner  frequency"  between  f  •  0.1  and  f  -  .5  Hz. 

The  patterns  show  the  effects  of  rupture  Propagation  at  frequencies  above^ 
the  corner  frequency*  due  to  contributions  fron  the  higher  order  P 
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indicated  by  equation  (16)  for  example.  Then  we  would  have  spectral  charac¬ 
teristics,  for  each  one  of  these  "events",  which  behave  as  previously  described, 
with  'the  complete  seismic  field  from  the  variable  stress  drop  event  being  a 
linear  superposition  of  all  the  "component  event"  spectra,  with  appropriate 
time  or  phase  delays  added  of  course. 

In  this  regard,  Figure  5  shows  a  model  of  the  San  Fernando  earthquake,  by 
Bache  and  Barker  (1977),  which  employs  a  multiple-event  superposition  of 
relaxation  theory  sources.  In  this  case  two  source  models  with  variable  stress 
drop  and  rupture  rate  were  employed.  The  near  field  model  predictions  for  the 
ground  motion  at  the  nearby  Pacoima  Dam  site  are  shown  to  be  in  good  first 
order  agreement  with  the  observations  at  this  location.  However,  in  view  of  the 
departure  of  the  predictions  from  the  observations  after  about  6  seconds  from 
the  start,  it  seems  evident  that  a  third  shallow  "event”,  intersecting  the  free 
surface,  could  be  added  to  the  model  to  give  a  more  complete  fit  to  the  entire 
observed  time  series.  Nevertheless,  the  first  six  seconds  of  the  observed  ground 
motion  is  predicted  reasonably  well  by  the  model. 

Ve  observe  that  the  detailed  nature  of  the  departure  of  the  predicted  and 
observed  ground  velocities  shows  that  the  true  velocity  is  more  complex  than 
that  predicted,  with  relatively  small  high  frequency  variations  in  the  recorded 
ground  velocity  that  are  only  fit  by  the  model  in  an  average  sense.  This  indi¬ 
cates  that  the  event  actually  has  fine  scale  variations  in  stress  drop  and  rup¬ 
ture  rate  that  involve  spatial  fluctuations  in  these  variables  on  a  scale  length  of 
hundreds  or  tens  of  meters,  as  well  as  on  the  scale  length  of  kilometers,  as 
represented  in  the  model. 

In  spite  of  the  lack  of  fine  scale  detail,  this  simple  "double  event”  model 
quite  accurately  predicts  the  observed  teleseismic  radiation  from  this  earth¬ 
quake.  Figure  6  shows  observed  far  field  seismograms  at  10  locations,  with  the 
theoretically  predicted  results  superimposed,  as  the  darker  lines,  on  the  obser¬ 
vations.  Comparison  of  the  upper  figure  results,  for  the  "event  1”  of  Figure  5 
alone,  with  those  in  the  lower  figure  where  both  "events”  are  included,  shows 
the  improvement  in  the  fit  to  the  data  when  the  two  event  model  is  used.  It  is 
also  quite  clear  that  a  third  "event”  at  shallow  depth  is  needed  to  obtain  a  good 
fit  to  the  data  beyond  10  to  15  seconds  after  the  first  energy  arrival.  In  any 
case  the  data  is  considered  to  be  well  fit  by  the  average  "double  event”  model 
for  the  first  10  to  15  seconds,  particularly  when  we  take  account  of  the  lack  of 
very  detailed  information  regarding  the  crustal  structure  at  the  various 
receiver  sites,  so  that  some  deviations  are  due  to  local  site  structure  inaccura¬ 
cies  rather  than  source  model  inaccuracies.  The  lower  frequency  content  and 
narrower  band  wjdtb  of  this  data,  compared  to  the  near  field  data  shown  in  Fig¬ 
ure  5,  allows  this  model,  which  clearly  does  not  account  for  fine  scale  variations 
in  stress  drop  and  rupture  rate,  to  provide  a  good  fit  to  the  far  field  data. 

A  somewhat  more  elaborate  model  made  up  of,  say,  three  such  "events” 
and  with  some  of  the  more  major  fine  scale  stress  drop-rupture  rate  fluctua¬ 
tions  included  could  undoubtedly  be  found  that  would  give  a  more  detailed  fit 
to  the  entire  wave  field  from  the  earthquake.  However,  for  the  purpose  of  this 
discussion,  it  is  sufficient  to  be  able  to  point  out  that  models  based  on  relaxa¬ 
tion  source  theory  give  rather  detailed  fits  to  observed  data  when  plausible 
earthquake  stress  drop  and  rupture  rate  variables  are  used,  and  secondly  that 
comparisons  with  well  recorded  near  and  far  field  data  show  conclusively  that 
both  large  scale  and  fine  scale  spatial  variations  in  the  basic  event  parameters 
(including  the  failure  zone  geometry)  occur  and  are  highly  significant,  particu¬ 
larly  for  the  high  frequency  radiation  fields  from  these  events.  The  example 
shown  bare  also  demonstrates  that  approximations  using  a  superposition  of 


Figure  5.  Veer  field  acceleration  from  the  San  Fernando  earthquake  recorded  at  the  Pacolna  Daa  site,  along  with 
the  corresponding  velocity  and  displacement  obtained  by  Integration.  The  dotted  lines  are  the  near  field 
predictions  obtained  from  the  "double  event"  relaxation  theory  source  model  Indicated  by  the  schematics  on 
the  right.  The  rupture  rate  and  stress  drop  at  the  failure  surface  mere  variable*  with  these  two  source 
variables  constrained  to  be  proportional  ' j  each  other.  (From  Bache  and  Barker,  1977) 
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71gur«  6.  Teleselsalc  signals  from  the  San  Farnando  Earthquake  (light  lines) 
compared  to  theoretically  predicted  signals  from  (1)  the  single  relaxation 
source  ("Event  1"),  near  the  hypocenter  (upper  figure)  and  (2)  the  multiple 
(double)  relaxation  soutce  nodcl  of  Figure  5.  All  predicted  and  observed  time 
sarlaa  are  scaled  to  unit  amplitude  for  the  first  cycle  of  the  direct  F  wave, 
with  the  "gain  ratio"  being  the  true  first  cycle  amplitude  ratio  of  synthetic 
to  observed,  (from  Bache  and  Barker,  1977) 
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relaxation  source  models  (each  with  an  approximate  solution  of  the  formal 
source  representation  theory  for  a  simplified  geometric  failure  zone)  can  be 
used  to  make  good  "first  order"  predictions.  (For  details  of  the  source  theory 
approximation,  as  well  as  the  wave  propagation  methods  used,  see  Bache  and 
Barker.  1977.) 

It  is,  however,  also  possible  to  compute  relaxation  source  theory  models 
with  spatially  heterogeneous  initial  stress  fields  specified.  Results  from  one 
such  computation,  by  Stevens  (19B1),  are  shown  in  Figure  7.  In  this  case  the 
rupture  zone  is  simply  specified  to  be  spherical  and  to  be  formed  "instantane¬ 
ously",  so  that  the  model  actually  conforms  to  an  explosion  produced  failure 
zone  in  a  heterogeneously  prestressed  medium.  Nevertheless,  the  features  of 
the  predicted  seismic  radiation  from  this  model,  that  are  solely  due  to  the 
departure  of  the  prestress  from  homogeneity,  are  easily  discerned.  In  particu¬ 
lar:  (9)  Both  the  P  and  S  wave  spectra  have  a  quadrupole  spectral  component 
that  is  the  same  as  that  for  a  source  of  the  same  geometry  but  with  uniform 
prestress  equal  to  the  spatial  average  of  the  non-uniform  prestress  field. 
Specifically,  this  quadrupole  field  component  has  a  flat  spectral  level  at  fre¬ 
quencies  below  the  comer  frequency  and  has  earner  frequencies  and  high  fre¬ 
quency  spectral  slopes  that  scale  in  the  same  manner  as  was  described  for  the 
uniform  prestress  case.  In  addition  however,  the  source  has  higher  order  mul¬ 
tipole  components  that  are  entirely  due  to  the  heterogeneity  of  the  initial  stress 
field  and  produce  interference  effects  which  are  superposed  on  the  quadrupole 
field  in  such  a  way  as  to  produce  peaked  spectra,  which  are  particularly 
apparent  near  nulls  in  the  quadrupole  radiation  pattern.  The  relative  strength 
of  the  higher  order  multipoles  is  directly  proportional  to  the  deviation  of  the 
heterogeneous  stress  level  from  the  mean  stress  and  the  contributions  to  the 
radiation  field  are  most  significant  at  and  above  the  quadrupole  spectrum 
comer  frequency  for  a  given  P  or  S  wave  type.  The  larger  the  spatial  extent  of 
the  deviation  in  the  initial  stress  from  the  mean,  the  lower  will  be  the  frequency 
at  which  there  will  be  a  significant  non-quadrupole  contribution  to  the  radia¬ 
tion  field. 

Figure  B  illustrates  how  a  strong  prestress  heterogeneity  can  produce  an 
"apparent  corner  frequency"  in  the  observed  spectra  from  a  tectonic  source 
which  is  more  representative  of  the  characteristic  dimension  of  the  zone  of 
high  stress  than  it  is  of  the  maximum  dimension  of  the  failure  zone.  This  exam¬ 
ple  is  again  from  Stevens  (19B1),  and  illustrates  the  radiation  induced  by  the 
creation  of  a  spherical  failure  zone  (e.g.,  as  might  be  produced  by  an  explosive 
shock  wave)  in  a  strongly  inhomogeneously  stressed  medium.  In  this  case  the 
prestress  is  created  by  insertion  of  a  fixed  ("locked")  dislocation  in  the  material 
at  a  variety  of  distances  from  the  spherical  failure  zone,  with  the  dislocation 
therefore  producing  the  initial  stress  field  and  the  relaxation  of  this  stress,  due 
to  creation  of  the  spherical  shatter  zone  with  vanishing  rigidity,  producing  the 
elastic  radiation.  As  is  shown,  when  the  dislocation  is  at  a  distance  of  ten  times 
the  shatter  zone  radius  (fio)  from  the  center  of  this  failure  zone,  then  the  pres¬ 
tress  in  the  vicinity  of  the  failure  region  is  fairly  uniform  and  a  "normal", 
predominantly  quadrupole  spectrum  is  produced  with  a  comer  frequency 
approximately  equal  to  the  rupture  rate  (here  taken  as  the  P  wave  velocity  in 
the  material)  divided  by  the  radius  of  the  failure  zone.  However,  when  the  dislo¬ 
cation  is  located  at  a  distance  of  1.1  /?0,  and  so  very  near  the  failure  zone  boun¬ 
dary,  then  the  quadrupole  spectrum  is  strongly  perturbed  by  the  addition  of 
higher  order  multipole  terms  arising  from  the  strong  spatial  variability  of  the 
initial  prestress.  These  higher  order  multipole  terms  are  seen,  at  least  in  this 
case,  to  produce  high  frequency  contributions,  at  and  above  the  quadrupole 
corner  frequency,  which  result  in  a  much  higher  apparent  corner  freuency. 
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This  "apparent"  corner  frequency  implies  a  characteristic  dimension,  when 
interpreted  in  terms  of  a  length  divided  by  rupture  rate,  which  is  of  the  order  of 
the  radius  of  the  high  stress  zone  surrounding  the  static  dislocation  rather 
than  the  (much  larger)  radius  of  the  failure  zone  itself.  Furthermore,  when 
both  the  spectral  amplitude  and  phase  are  used  to  obtain  the  time  domain 
seismograms  corresponding  to  this  situation,  it  is  found  that  a  large  high  fre¬ 
quency  pulse  appears  with  a  time  of  arrival,  determined  by  the  phase  spectrum 
and  associated  groups  delay,  that  corresponds  to  a  source  of  energy  in  the 
vicinity  of  the  dislocation.  Therefore,  even  though  the  dislocation  displacement 
discontinuity  is  fixed  for  all  time,  the  energy  is  clearly  found  to  arrive  from  the 
dislocation.  This,  on  reflection,  would  be  expected  from  a  volume  relaxation 
source  where  the  energy  is  released  from  the  medium  itself,  and  would  be  most 
evident  when  strong  heterogeneities  in  the  prestress  exist. 

fVom  the  standpoint  of  predicted  characteristics  of  the  radiation  from 
earthquakes,  it  is  therefore  important  to  keep  in  mind  that:  (10)  Both  the  Pand 
S  waves  can  have  spectral  properties  that  reflect  the  characteristics  of  strong , 
nearby,  stress  concentrations,  in  particular  very  high  stress  levels  can  produce 
large  energy  arrivals  within  the  P  and  S  wave  trains  which  produce  apparent 
corner  frequencies  of  the  composite  spectra  of  the  train  that  are  more  charac¬ 
teristic  of  the  dimension  of  the  stress  concentration  than  of  the  failure  sane 
dimension.  However,  the  energy  arrival  (s)  producing  this  'perturbation  "  in  the 
Par  S leave  amplitude  spectra  have  travel  times  appropriate  to  energy  released 
from  the  location  of  the  stress  concentration,  which  may  not  be  the  same  as  the 
location  of  the  point  of  initiation  of  failure  (the  hypocenterj,  and  so  such  energy 
may  frequently  arrive  in  the  P  or  S  wave  coda  after  the  first  arriving  energy, 
although  it  could,  in  fact,  constitute  the  first  arriving  pulse  if  the  hypocentral 
■  sane  was  Initially  an  intensely  stressed  sons.  ■  - 

Scaling  laws  for  Earthquakes  Based  on  Theoretical  Results 

The  predicted  properties  of  the  seismic  radiation  from  these  dynamical 
models  can  also  be  obtained  analytically  if  we  consider  the  asymptotic  behavior 
of  relaxation  source  model  radiation,  such  as  described  by  equations  (34) 
through  (37).  This  approach  provides  quantitative  relations  which,  while 
approximate,  are  valid  to  first  order  and  can  be  used  as  scaling  laws  for  the 
radiated  spectrum  for  earthquakes.  That  is,  we  can  scale  the  spectra  for  one 
particular  tectonic  source  into  the  spectra  for  a  different  source  with  a  dif¬ 
ferent  stress  drop,  rupture  length  and  rupture  velocity  using  these  relations. 

Since  the  theory  predicts  results  that  are  different  for  the  compressional 
and  shear  waves  radiated  by  tectonic  sources,  we  will  treat  the  spectra  for 
these  wave  types  separately.  Most  of  the  following  results  are  from  Archambeau 
(1968,  1972)  and  Minster  (1978),  and  details  omitted  here  can  be  found  in  these 
references,  particularly  in  Minster  (1978). 

In  order  to  simplify  the  results,  without  losing  any  essential  generality,  we 
consider  the  case  in  which  the  prestress  field  is  homogeneous  and  pure  shear, 
with  only  the  initial  stress  component  nonzero.  Then,  using  the  expression 
(34)  for  the  radiation  in  the  limit  of  low  frequencies  and  in  the  far  field  (so 
kpr  »  1),  we  have  for  the  compressional  (P)  wave  field  (Minster,  1978): 

Urn  HP(r.  u)  ~  J^77-2|o)  ^7  ,V  Bin  2*  cos*  <37> 

where  the  asymptotic  result  has  been  expressed  in  terms  of  the  radial  com¬ 
ponent  of  the  displacement  spectrum  Here  a  is  the  Poisson's  ratio  for  the 
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material  in  the  region  eurrounding  the  failure  cone,  ft  ie  the  elastic  rigidity  in 
this  zone,  9$  the  prestress,  vp  the  compression  wave  velocity  et  the  source. 
kp  m  o/vp  is  the  wave  number  and  L  is  the  maximum  dimension  (or  length")  of 
the  failure  sone.  This  limiting  form  of  the  low  frequency  far  field  radiation 
•hows  that  the  field  has  a  simple  quadrupole  form  which  is  directly  proportional 
to  the  prestress  and  the  cube  of  the  rupture  length. 

Using  this  low  frequency  result  along  with  asymptotic  results  valid  at  high 
frequencies,  we  can  define  a  "P  wave  corner  frequency"  (or  olp*  =  Znf^) 
ms  the  point  where  the  two  asymptotic  spectral  results  are  equal.  In  addition, 
we  can  epecify  the  high  frequency  behavior  of  the  spectrum  from  the  high  fre¬ 
quency  asymptotic  results.  Knee  the  high  frequency  asymptotic  relations  are 
rather  complex  we  simply  summarize  their  consequences  when  used,  along  with 
(37).  to  define  a  "corner  frequency"  and  the  high  frequency  behavior  of  the  P- 
wave  spectrum,  at  frequencies  above  Ve  have: 

* — hsrj  (3B) 

for  the  angular  "corner  frequency"  relation,  and 


lim  «*•  O  -tt  ;  when  Up  <  vP 

h»i  tr 


lim  ~  0  -V  ; 


when  Up  fc  vp 


for  the  very  high  frequency  asymptotic  behavior  of  the  far  field  P  wave  spec¬ 
trum  as  a  function  of  frequency.  Here  Up  is  the  rupture  rate  and  the  notation 

O  Jjjj.  for  example,  is  to  be  read  as  "the  order  of  1/  u*".  The  behavior  at  lower 

frequencies,  but  such  that  o  Jfc  oip\  where  the  frequency  is  near  but  somewhat 
larger  than  the  corner  frequency,  is  such  that: 

lim  ~  0  I -V| ,  for  Up  <  vj»  and  Cfe  fc  vp  (40) 


Thus,  the  P  wave  spectrum  decreases  with  a  slope  near  o'*  at  frequencies  above 
t»lPi  and,  when  the  rupture  rate  Up  is  less  than  the  compressions!  velocity,  then 
this  slope  increases  to  o'*  at  yet  higher  frequencies.  As  the  rupture  rate 
approaches  the  congressional  velocity  in  value,  the  frequency  at  which  the 
spectrum  assumes  a  o~*  behavior  moves  to  higher  frequencies,  such  that  when 
the  rupture  velocity  actually  reaches  or  surpasses  the  compressional  velocity 
»  (*.y„  as  for  failure  caused  by  a  shock  wave),  then  the  frequency  at  which  o'* 
behavior  occurs  is  infinite.  Further,  the  spectral  decay  beyond  the  corner  fre¬ 
quency  ojft  is  strongly  influenced  by  the  ratio  of  Up  to  vp,  and  is  essentially 
independent  of  the  rupture  dimensions.  When  the  rupture  rate  is  much  less 
than  the  shear  velocity,  and  therefore  very  much  less  than  the  compressional 
velocity,  then  the  spectral  decay  ranges  from  «”*  to  o~*  over  a  broad  frequency 
band  and  the  corner  frequency  predicted  by  equation  (38)  is  not  very  accurate. 
This  behavior  is  illustrated  in  Figure  3,  and  was  also  noted  earlier. 

These  results  are  appropriate  for  a  uniform  prestress  field  and,  as  was 
illustrated  and  discussed  earlier,  prestress  inhomogeneities,  particularly  strong 
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stress  concentre tiona,  can  affect  the  apectrum  near  and  above  the  corner  fre¬ 
quency  so  aa  to  cause  peaking  and  corner  frequency  shifts.  These  effects  will 
nearly  always  arise  from  later  arriving  energy  within  the  P  wave  train,  and  so 
the  first  arriving  P  wave  pulse  will  typically  have  a  simpler  spectrum,  with 
characteristics  most  closely  matched  by  the  results  predicted  in  (37)-(39).  It  is 
only  when  the  spectrum  of  the  entire  P  wave  train  is  considered  that  compli¬ 
cated  Interference  of  arrivals  in  the  train  produces  strong  perturbations  in  the 
total  P>wave  spectra. 

Using  these  results  we  can  construct  approximate  scaling  laws  for  the  P 
wave  spectra  between  different  sources  having  different  fault  length,  stress 
drop  and  rupture  rate  parameters.  These  relations  are  accurate  and  useful 
when  the  rupture  rate  is  less  than,  but  fairly  near,  the  shear  velocity  in  the 
medium,  which  appears  to  be  the  usual  case  for  failure  in  the  earth,  as  was 
noted  earlier.  The  relation  (37)  also  holds  for  the  "supersonic"  failure  rate, 
when  Ur  at  vp.  However,  the  corner  frequency  estimate  is  less  accurate  for  a 
very  low  rupture  rate  process,  and  the  usefulness  of  the  relation  in  (3B)  for  the 
corner  frequency  is  minimal  in  this  case.  Nevertheless,  with  these  provisions, 
the  results  given  in  (37)  through  (39)  provide  the  following  scaling  relations  for 
the-P  wave  spectra  from  relaxation  models  of  earthquakes: 

(1)  The  spectral  level  of  the  for  field  P  wave  at  low  frequencies,  below  the 
comer  frequency,  is  independent  of  frequency  and  rupture  rate  and  this 
flat  level  scales  with  fault  dimension  according  to: 

(41) 

where  A^  («)  denotes  the  low  frequency  spectral  amplitude  level  of  a 
reference  event,  while  A^{u)  is  the  flat  spectral  level  of  an  event  with  a 
different  failure  zone  dimension,  but  with  oil  other  source  parameters  the 
earns  as  the  reference  event.  This  result  is  most  accurately  obeyed  when 
the  failure  zone  has  a  second  dimension,  W,  comparable  in  size  to  L,  that  is 
of  the  same  order  of  magnitude.  Further,  it  is  assumed  that  W  changes 
proportionally  with  L.  If  this  is  not  the  case,  then  a  somewhat  more  accu¬ 
rate  scaling  relation  is: 

A^Xu)  f  WL*  ) 

In  the  case  of  a  spherical  failure  zone,  of  radius  R,  then 
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The  spectral  level  of  the  entire  P  wave  spectrum  scales  directly  with  the 
prestress,  and  so: 


|S>)L  ,  L»»>L 

mP’wi  tin 


(42) 


where  |u^|  denotes  the  entire  P  wave  amplitude  spectrum  for  a  reference 
event  while,  |u^|  Is  the  spectrum  for  an  event  having  a  different  homo¬ 
geneous  (shear)  prestress  magnitude,  but  with  all  other  source  parameters 
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nks  ram*  os  tho so  for  tho  reference  event,  (Here  it  is  assumed  that  the 
shear  stress  drop  at  the  failure  boundary  is  total  so  that  the  scaling  rela¬ 
tion  is  given  in  terms  of  the  initial  prestress.  If  the  stress  drop  at  the 
boundary  of  the  failure  tone  is  partial,  then  the  stress  quantities  in  (42) 
denote  shear  stress  changes  across  the  failure  sone  boundary.) 

The  corner  frequency.  fP»  for  the  compressions!  wave  spectrum  from  the 
source  scales  directly  os  the  ratio  of  rupture  velocity  to  failure  sons  length, 
and  in  particular  for  rupture  rates  near  shear  velocity  in  the  medium  (the 
most  common  situation)  then: 


»mhrj  *882|t-|; 


£*5  and  o-  — 


Further,  event  spectrum  scaling  relative  to  a 
"corner  frequency”  denoted  by  Wfjft,  is  given  by 
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reference  source,  with 


Ve  note  that  once  the  scaling  laws  in  (41)-(43)  have  been  appropriately 
applied  to  the  spectrum  of  a  reference  event,  then  the  proper  flat  spectral  level 
at  low  frequencies  will  have  been  established  for  the  new  event,  along  with  its 
appropriate  corner  frequency.  Then,  for  P-wave  spectra,  there  will  be  a  rather 
short  frequency  band  above  where  the  spectrum  decays,  roughly  as  w-*, 
before  achieving  the  predicted  u~*  decay  at  higher  frequencies.  However,  since 
the  rupture  velocities  of  natural  events  is  such  that  Up*vs<  vp,  so  that  Up  is 
significantly  less  than  Vp.  then  the  spectrum  achieves  a  slope  quite  rapidly 
with  increasing  frequency.  (This  behavior  is  illustrated  in  Figure  3  by  the  P- 
wave  spectra  with  rupture  velocity  Up  =  3.45  km/sec  and  Up  =  3.0  km/sec.) 
Thus,  assumption  of  a  ratio  of  the  rupture  rate  to  compressional  wave  velocity 
establishes  the  spectral  decay  for  frequencies  above  the  corner  frequency  and 
since  this  ratio  is  always  significantly  less  than  unity  for  natural  earthquake  (or 
spontaneous)  failure,  then  the  first  equation  in  (39)  applies  and  the  u~9  spec¬ 
tral  decay  holds  over  nearly  all  of  the  frequency  band  where  /  >  fp,  In  this 
way  then,  the  first  order  spectrum  of  the  compressional  wave  field  can  be  esta¬ 
blished  with  reasonable  accuracy  over  the  entire  frequency  range. 

As  noted  earlier  however,  when  the  rupture  rate  is  controlled  by  a  high 
speed  shock  wave,  such  as  is  produced  by  an  explosion,  then  the  rupture  rate  is 
close  to  the  compressional  wave  velocity  and  a  «“*  spectral  decay  rate  beyond 
the  corner  frequency  will  occur. 

Therefore,  the  scaling  laws,  plus  knowledge  or  assumption  of  the  rupture 
rate,  can  be  used  to  provide  an  estimate  of  the  entire  P-wave  spectrum  gen¬ 
erated  by  the  relaxation  of  tectonic  stress  in  the  medium  surrounding  a  natural 
or  explosion  induced  failure  cone.  Of  course,  in  view  of  the  previously  discussed 
predictions  and  observations  concerning  prestress  concentration  effects,  the 
spectrum  produced  by  scaling  will  not  include  the  perturbative  effects  of  inho¬ 
mogeneous  prestress.  However,  the  scaled  spectrum  should  be  a  good  average 
(or  smoothed)  representation  in  general,  and  only  in  (rare)  cases  of  quite 
strong  prestress  concentrations  will  significant  deviations  occur.  Further,  even 
when  strong  prestress  concentrations  do  occur,  they  will  most  often  be  mani¬ 
fested  in  the  later  arriving  energy  within  the  P  wave  train  and  will  not  signifi¬ 
cantly  affect  the  first  arriving  energy.  Thus,  for  predictions  of  the  first  arrival 
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energy,  the  scaled  spectra  should  usually  be  quite  accurate. 

The  structure  and  form  of  the  S  wave  spectrum  produced  by  a  tectonic 
source  has  some  important  differences  from  the  P-wave  spectrum.  In  particu¬ 
lar,  the  expression  comparable  to  the  far  field  P  wave  spectrum  at  low  frequen¬ 
cies,  given  by  equation  (37),  is: 

Um  ai‘J(r.U)  ~  [£] •“*  <44> 


Here  represents  either  the  d  or  p  displacement  component  and  vs  is  the 
shear  velocity.  This  field  has  a  quadrupole  form  and  also  has  a  flat  spectral 
amplitude,  as  did  the  P  wave.  However,  the  ratio  of  the  iou  frequency  ampli¬ 
tude  of  the  S  wave,  A^(«),  to  the  same  amplitude  function  for  the  P  wave  is 
easily  seen,  from  (37)  and  (44),  to  be  given  by: 


■A^(cj)  _  Vp_  1  —  a  _  f  vp  j3 
A^(«)  “  vs  1  -  2a  “  [  vs  j 


(45) 


Therefore,  a  particular  tectonic  source  will  have  a  low  frequency  S  wave  ampli¬ 
tude  spectrum  that  is  larger  than  the  P  wave  spectrum  by  a  factor  equal  to  the 
cube  of  the  ratio  of  the  P  to  S  wave  velocities  in  the  medium.  For  a  typical  solid 
this  factor  is  about  3V3  and  so  represents  a  very  significant  difference.  (A 
difference  in  P  and  S  spectra  of  about  this  order  of  magnitude  is  clearly 
observed  for  earthquakes,  so  that  this  prediction  is  well  supported  by  direct 
observations.) 


In  addition,  the  predicted  S-wave  spectra  can  be  shown  to  have  a  corner 
frequency  given  (approximately)  by: 


2  Ur 
L 


(46) 


Thus  the  corner  frequency  for  the  S-wave  spectra  is  different  than  that  for  P 
waves,  and  is  lower  because  of  the  dependence  on  the  ratio  of  vs  to  Ur.  as  com¬ 
pared  to  the  P  wave  corner  frequency,  which  has  the  same  form  but  depends  on 
the  (larger)  raUo  of  to  Ur.  In  the  usual  case,  when  the  rupture  rate  is  near 
the  shear  velocity  value  in  the  material,  then 


Ur  &vs 


this  value  is  about  30%  lower  than  the  P  wave  corner  frequency  for  the  same 
event. 


The  high  frequency  asymptotic  behavior  of  the  S  wave  spectra  is  given  by: 

(47) 


Km  ~  O  j~j  ;  when  Ur  <  vs 

Urn  |ff$|  ~  0  |~J  ;  when  Ur  *vs 


Because  of  the  fact  that  the  asymptotic  behavior  is  different  for  the  different 
ranges  of  rupture  velocity,  and  since  the  rupture  rate  is  usually  of  the  order  of 
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the  shew  vm  Telocity  (it.,  Efe  i*  only  eligbtly  leu  or  equal  to  the  shear  Telo¬ 
city).  then  the  S  wave  spectra  trill  have  a  very  broad  frequency  range  within 
which  the  spectral  decay  caries  approximately  as  jT*.  This  is  in  contrast  to  the 
P  wave  spectra,  for  which  the  spectral  band  within  which  the  o’*  decay  applies 
is  Just  above  of  and  is  quite  narrow,  while  over  most  of  the  range,  where 
«i  >  ojfK  the  decay  Is  o~*.  Thus,  when  the  rupture  rate  is  close  to  the  shear 
velocity  In  the  medium,  then  the  S  wave  decay  is  cT*  out  to  quite  high  frequen¬ 
cies.  relative  to  of^,  before  it  assumes  a  o~*  decay.  When  the  rupture  rate  is 
larger  than  the  S  wave  velocity,  then  the  S  wave  spectral  decay  is  «“*  for  all  fre¬ 
quencies  above  Observationally,  the  S  wave  spectrum  should  appear  to 

have  a  o~*  decay,  while  the  P  wave  spectrum  will  decay  as  o’*  in  essentially  all 
eases  except  for  very  slow  rupture  rates,  or  for  rupture  rates  near  or  at  the 
compression^  velocity.  The  latter  situation  probably  only  applying  for  explo¬ 
sive  shock  induced  failure. 


Based  on  the  results  in  (44)-(47)  the  following  scaling  laws  generally  apply 
for  S  wave  spectra: 

(1)  The  spectral  level  of  the  far  field  S  wave  at  law  frequencies,  below  the 
comer  frequency,  is  independent  of  the  rupture  rate  and  frequency,  and 
the  flat  level  scales  with  fault  dimension  according  to: 


(4B) 


where  A^iu)  denotes  the  low  frequency  spectral  level  of  a  reference 
event,  while  A^(a)  is  the  low  frequency  spectral  level  of  an  event  with  a 
different  failure  zone  dimension,  but  with  all  other  source  parameters  the 
same  as  the  reference  event.  When  the  second  largest  failure  zone  dimen¬ 
sion.  W.  is  not  of  the  same  order  as  L.  or  does  not  change  proportionally 
with  L  tor  the  events,  then  a  more  accurate  relation  is: 

A^(cj)  [  WL*  \ 

aFHo) 


.For  spherical  rupture  zones,  then  (48)  holds  with  L  the  diameter  of  the 
spherical  region.  Finally,  the  low  frequency  S  wave  spectral  level  is  related 
to  the  P  wave  spectral  level  by 


where  vp  and  v$  are  the  compressional  and  shear  velocities  in  the  medium 
surrounding  the  failure  zone. 

The  spectral  level  of  the  entire  5  wavs  spectrum  scales  directly  with  the 
prestress,  and  so: 


where  6$^  denotes  the  entire  S  wave  amplitude  spectrum  for  a  reference 
event,  while  ufS)  is  the  spectrum  for  an  event  having  a  different  homogene¬ 
ous  (shear)  prestress  magnitude,  or  a  different  homogeneous  stress  drop, 
but  with  all  other  source  parameters  Vie  same  as  those  for  the  reference 
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(3)  The  corner  frequency,  /j*1,  for  the  sheer  wave  spectrum  from  the  source 
•cafes  diroctly  as  tho  ratio  of  rupturo  voloctty  to  faCuro  so no  i anyth.  Par 
rupturo  raise  near  tho  ohoor  velocity,  as  is  usually  tho  e ass; 

fF’n  ix]  wr  -  •4*»  |t-]:  ■ 1/4 


The  ratio  of  P  to  S  wave 


corner  frequencies  for  a  particular  event  is: 

fP  Mf* 

TF’RJ 


Further,  event  spectrum  scaling  relative  to  a  reference  source,  with 

"corner  frequency"  denoted  by  Mffr,  is  given  by 

wk-lwM  <50) 

Comments  concerning  the  effects  of  strongly  inhomogeneous  prestress  (or 
stress  drop)  on  P  wave  spectra  apply  with  equal  force  for  the  S  wave  spectra  as 
well.  However,  in  both  cases  the  scaling  laws  and  asymptotic  relations  will  ordi¬ 
narily  hold  in  an  average  sense,  so  that  the  relations  given  here  are  quite  use¬ 
ful. 

Figure  9  provides  a  graphical  summary  of  the  predicted  spectral  charac¬ 
teristics  for  the  eompressional  and  shear  waves  from  earthquakes.  In  this  fig¬ 
ure  the  prestress  (or  stress  drop),  denoted  as  |e|,  and  the  failure  length.  L,  are 
left  arbitrary,  so  that  the  ebsolute  amplitude  scale  is  in  arbitrary  units.  How¬ 
ever.  the  rupture  rate  is  essumed  to  be  close  to  the  shear  velocity  in  the 
medium,  which  is  the  most  commonly  occurring  situation  in  the  earth,  while  the 
ratio  of  fault  length  to  rupture  velocity  is  taken  to  be  unity.  (Thus,  for  exam¬ 
ple.  since  Cfe  at  «•  and  with  115  at  3.B  km/sec,  then  L  m  3.B  km.)  Finally,  the 
Poisson's  ratio  tor  the  material  is  taken  to  be  1/4,  which  is  also  typical.  The  low 
frequency  amplitude  asymptote  for  the  far  field  P  wave  has  been  normalized  to 
unity  in  the  plot. 

With  the  choice  vs  *=  3.B  km/sec.  then  vp  *  6.4  km/sec,  Ur  -  3.6  km/sec 
and  L  *  3.B  km  are  required,  under  the  assumptions  used,  to  obtain  the  spectra 
in  the  figure.  Using  the  relation  shown  for  Ap(o),  with  a  prestress  of  130  bars 
(1.3  x  10*  dynes/cm*)  and  a  rigidity  of  3  x  1011  dyne /cm8,  so  that 
|o|//a  =  4.3  x  10-4  is  the  strain  magnitude,  we  have: 

Ap{u)  m  10*  cm  sec 

Since  the  P  wave  displacement  spectrum  at  a  radial  distance  |r|  is  given  by 

■  J V  I  ^p(w)  n  . _ 


|  «/!(„)  | 


sin  2d  cosy 


from  equation  (37),  then  the  low  frequency  spectral  level  at  a  distance  of  r 
km  from  the  3.B  km  failure  sone  is 

IttPt—  *  1  cm -esc;  r  *  1 0km 


. . . .  .  I  J  J  I  ■ 

—  '  vz  l 


•08 


Hence  if  the  amplitude  scale  in  the  figure  is  in  cm-sec,  then  the  spectra  shown 
are  appropriate  at  a  distance  of  10  km  from  a  3.8  km  fault,  with  failure  occur¬ 
ring  at  a  uniform  rupture  rate  of  3.8  lcm/sec  in  a  medium  prestressed  to  130 
bars.  The  elastic  velocities  vs  =  3.8  km/sec,  vp  =  8.4  km/sec  and  rigidity 
ft «  3  x  1011  dynes/cm*  are  appropriate  to  the  middle  crust  of  the  earth  in  most 
areas  (*.p..  from  10  to  about  25  km  depth).  Further,  with  the  rupture  rate  near 
the  shear  velocity  and  the  ratio  of  rupture  rate  to  fault  length  unity,  then  the 
corner  frequencies  of  the  spectra  will  be  at  .662  Hz  for  the  P  wave  spectra  and 
at  .459  Hz  for  the  S  wave  spectra.  (However,  these  corner  frequency  values  do 
not  depend  on  the  absolute  values  of  L,  Ur  and  vs,  but  only  on  their  ratios.) 

figures  10a.  and  10b.  illustrate  an  application  of  the  scaling  laws  to  P  and 
8-wave  spectra,  for  the  case  in  which  only  the  rupture  length  varies.  Only  the 
far  field  spectra  are  indicated  and  only  the  smoothed  form  of  the  spectra  are 
indicated.  Also  it  is  assumed  that  the  prestress  field  is  nearly  homogeneous, 
with  a  relatively  small  stress  concentration  near  the  initial  rupture  point  giving 
the  slightly  peaked  spectra  Indicated.  (The  "peaking"  in  the  spectra  will  be  a 
function  of  azimuth  in  general  and  only  a  typical  azimuth  is  indicated  here.)  In 
any  case,  for  rather  "small"  rupture  lengths  with  L  £  10  km.  then  we  expect  the 
rupture  surface  areas  to  be  proportional  to  L2  and  that  the  L9  scaling  law  for 
amplitudes  can  be  applied.  Since  the  corner  frequency  for  both  P  and  S  waves 
scales  inversely  with  L.  then  as  a  consequence  of  these  two  scaling  laws,  the 
locus  of  corner  frequency  points  will  lie  along  a  straight  line  with  a  f~9  slope  in 
these  plots. .  Since  the  P  wave  spectrum  has  a  /"*  high  frequency  asymptote  for 
rupture  rates  near  the  shear  velocity,  then  the  spectra  for  P  waves  will  also 
approach  and  lie  along  this  "scaling  line"  at  frequencies  above  the  corner  fre¬ 
quency.  Therefore,  all  the  events  will  have  the  same  P  wave  spectral  levels  at 
high  frequency,  since  they  all  approach  the  /“•  scaling  line.  However,  S  wave 
spectra  have  a  /"*  asymptotic  behavior  over  a  wide  frequency  range  above  the 
corner  frequency,  and  so  the  S-wave  spectra  will  not  approach  the  scaling  line 
except  at  extremely  high  frequencies,  which  are  usually  outside  the  range  of 
seismic  observation. 

Figure  11  illustrates  the  applicable  spectra  scaling  for  variable  average 
stress  drops  or  changes  in  the  average  prestress  level.  Here  we  assume  that 
the  prestress  (or  stress  drop)  is  quite  homogeneous  and  that  the  changes  in 
stress  are  confined  to  changes  in  the  average  value.  For  such  changes,  both 
the  entire  P  and  S  wave  spectra  scale  directly  with  the  magnitude  of  the  stress 
and  hence  the  scaling  line  of  the  corner  frequencies  is  vertical.  Here  we  show 
only  the  far  Held  P  wave  spectra  since  the  S-wave  spectral  scaling  is  identical. 
In  addition,  the  near  Held  spectral  contributions  also  scale  directly  with  the 
stress  magnitude  and  so  behave  exactly  like  the  far  field. 

The  vertical  scaling  line  for  prestress  or  stress  drop  changes  will  be  modi¬ 
fied  if  the  rupture  rate  is  proportional  to  the  stress  drop  or  prestress.  This 
possibility  was  discussed  earlier  and  is  quite  likely  for  failure  processes  in  the 
earth.  If  this  is  the  case,  then  the  stress  scaling  must  be  combined  with  the 
rupture  velocity  scaling  law  illustrated  in  Figure  12.  As  shown  in  Figure  12.  the 
scaling  line  for  the  comer  frequencies  for  P  wave  spectra  is  a  horizontal  line. 
Further,  the  high  frequency  spectral  slope  can  change  from /"•,  for  Vp  £  vs.  to 
/-*,  for  Up  at  vp.  Also  for  very  low  rupture  rates  there  is  a  frequency  range 
above  the  corner  frequency  within  which  the  spectral  slope  is  low,  varying  from 
/"*  to  /"•  in  a  continuous  fashion. 

If  we  ware  to  combine  the  rupture  rate  scaling  with  the  stress  drop  scaling, 
as  would  be  required  for  rupture  rates  proportional  to  stress  drop,  then  we 
would  expect  to  see  a  corner  frequency  locus  defining  a  scaling  line  that  had  a 


Pi  gu re, 10  a .  Illustration  of  an  application  of  the  scaling  laws  to  P  wave  spectra  for  variable 
rupture  length  (L),  with  all  other  source  parameters  fixed.  The  corner  frequency 
an  f~3  locus,  defining  a  "Scaling  Line".  Scaling  assumes  rupture  surface  area 
increases  as  ("Small"  L  case) . 


FREQUENCY,  f 

Illustration  of  an  application  of  the  scaling  laws  to  S-vave  spectra  for  variable 
rupture  length  (L) ,  with  all  other  source  parameters  fixed.  The  corner  frequency 
traces  an  f-3  locus,  defining  a  "Scaling  Line"..  Scaling  assumes ■ rupture  surface 
area  increases  as  L*.  ("Small"  L  case) 


Figure  11.  Illustration  of  an  application  of  the  scaling  laws  to  P-wave  spectra  for  variable 
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Figure  12.  Application  of  scaling  laws  to  P  wave  spectra  asymptote  for  variable  rupture  velocity  (U_),  with 
all  other  source  parameters  fixed.  The  locus  of  corner  frequencies  defines  a  horizontal  ’'Scaling 
Line".  Bote  that  the  high  frequency  asymptote  for  Ur  ~  vp  is  f”2,  while  for  lower  Ur  values 


positive  slope  /*  or  /*,  say)  rather  than  a  vertical  scaling  line,  and  the 
spectra  for  higher  stress  drop  events  would  have  a  higher  corner  frequency 
than  those  for  lower  stress  drops.  Further,  one  could  expect  to  observe  slight 
ehahges  in  high  frequency  spectral  slopes.  As  noted,  this  behavior  is  more 
likely  than  not  to  occur  in  the  earth,  and  should  certainly  be  considered  when 
applying  the  scaling  relations. 

The  scaling  relations  in  Figure  12,  as  well  as  in  Figure  11,  are  applied  only 
to  far  field  P  wave  spectra,  however,  the  situation  is  nearly  Identical  for  the  S- 
wave  spectra  except  for  the  fact  that  the  high  frequency  asymptotic  behavior  is 
different  for  the  S-wave.  This  difference  does  not  lead  to  any  fundamental 
changes  in  the  S-wave  spectral  scaling  results  however. 

The  applications  of  the  scaling  laws  in  Figures  9-12  have  been  limited  to 
rather  "small"  tectonic  sources,  that  is  to  events  with  failure  tone  dimensions 
that  are  small  enough  that  we  can  expect  that  the  surface  area  of  the  failure 
sone  will  scale  as  L *.  This  probably  applies  to  events  with  L  £  10  km.  However, 
for  much  larger  events,  we  certainly  do  not  expect  that  the  area  is  proportional 
to  the  square  of  the  length  or  maximum  dimension  of  the  failure  zone.  Thus  for 
large  events,  with  rupture  zones  of  the  order  of  hundreds  or  thousands  of 
kilometers,  the  second  largest  dimension,  W,  would  be  expected  to  be  nearly 
fixed  in  value  and  independent  of  length  of  these  large  failure  zones.  In  this 
case,  the  surface  area  of  the  failure  zone  would  scale  as  L.  rather  than  as  £*, 
and  the  scaling  relation  for  the  low  frequency  amplitude  would  be  proportional 
to  L*  rather  than  to  L9.  This  would  have  the  effect  of  producing  a  scaling  line, 
for  the  corner  frequencies  of  both  P  and  S  wave  spectra,  that  varied  as  f  with 
failure  zone  length  for  the  very  large  events.  Thus,  while  this  line  varies  as  / 
for  small  events  with  L  £  10  km,  it  would  be  expected  to  vary  as  /  "*  for  very 
large  events,  certainly  for  events  with  L  Jfc  100  km.  In  the  intermediate  range, 
where  10  £  L  £  100  km.  the  variation  of  the  slope  of  this  scaling  line  should  be 
between  f~9  and  /“*,  with  the  variation  being  continuous  between  the  two, 
reaching  these  slopes  smoothly  at  the  extreme  values  of  L.  near  100  km  and  10 
km,  respectively. 

A  consequence  of  this  scaling  of  the  spectra  with  failure  zone  dimension  is 
that  the  P  and  S  wave  spectral  levels  will  cease  to  increase,  at  frequencies 
above  the  corner  freuency,  even  though  the  failure  dimension  L  continues  to 
Increase.  This  occurs,  as  already  shown  in  Figure  10a,  for  P  wave  spectra  at 
even  the  smallest  L  values.  On  the  other  hand,  as  shown  in  Figure  10b,  the  S 
wave  continues  to  have  larger  high  frequency  spectral  levels  with  increasing  L, 
for  the  lower  L  value  range.  However,  for  L  values  that  are  large,  the  amplitude 
scaling  law  varies  as  £*,  rather  than  as  L9,  and  the  scaling  line  has  a  /  ~*  slope, 
so  that  in  this  large  L  range  the  high  frequency  S  wave  spectra  will  cease  to 
increase  with  increasing  L.  That  is,  as  in  the  case  for  the  P  wave  spectra  in  the 
low  L  range,  there  is  "spectral  saturation"  in  the  high  frequency  range  for  S 
wave  radiation  and  the  high  frequency  S  wave  spectral  level  will  not  increase 
with  increasing  L. 

Am  will  be  emphasized  later,  this  high  frequency  "saturation"  for  P  waves 
leads  to  a  out-off  in  the  body  wave  magnitude  (m*)  measured  from  P  waves, 
where  increases  in  the  rupture  dimension  L  do  not  result  in  an  increase  in  the 
body  wave  magnitude  measured  at  1  Hz.  This,  of  course,  arises  from  the  fact 
that  the  body  wave  magnitude  is  defined  as  the  logarithm  of  the  one  second 
period  P  wave  amplitude.  Thus,  when  the  corner  frequency  is  significantly  less 
than  1  Hz,  then  the  1  Hz  point  will  be  on  the  /“•  asymptote  of  the  P  wave  spec¬ 
trum  and  all  events  with  larger  failure  zone  lengths  will  still  have  (very  nearly) 
the  same  spectral  amplitude  at  1  Hz.  Hence  there  will  be  no  change  in  the 
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magnitude  beyond  this  limiting  value  and  one  will  observe  a  body  wave  magni¬ 
tude  cut-off.  Observationally,  this  clearly  occurs  for  body  wave  magnitudes 
somewhat  above  magnitude  6.  (This  corresponds,  roughly,  to  a  failure  zone 
dimension  L  of  the  order  of  10  to  15  km)  On  the  other  hand,  when  the  failure 
cone  dimensions  become  very  large,  with  L  of  the  order  of  100  km  and  more, 
then  both  P  and  S  wave  high  frequency  spectra  saturate  above  their  comer  fre¬ 
quencies.  In  this  case  measurements  of  a  surface  wave  magnitude  (Us)  defined 
In  terms  of  the  logarithm  of  the  (vertical  component)  Rayleigh  wave  amplitude 
at  .05  Hz,  will  also  cease  to  increase  with  increasing  failure  zone  dimension. 
Thus,  there  is  a  surface  wave  magnitude  cut-off  which  will  apply  to  very  large 
tectonic  events.  Such  a  cut-off  has  not  been  conclusively  verified  observation- 
ally.  but  should  appear  at  an  Us  value  somewhat  above  Us  *  B.  For  such  large 
events  then,  neither  m*  nor  Us  will  give  a  measure  of  the  size  of  the  event,  or 
its  energy. 


